A maximum exact match parallel lookup algorithm for accelerating third generation RNA sequencing data alignment

Abstract

The application discloses a maximum exact match parallel search algorithm for accelerating third-generation RNA sequencing data alignment, and is suitable for accelerating an alignment algorithm designed on the basis of an "index-seed-link" idea. The algorithm firstly performs BWT transformation on reference sequences of a genome to construct a corresponding FM index; then, all prefixes of the reference sequences of the genome are calculated, and the prefixes are sorted, so that a longest common prefix array is created and sampling is performed. In order to accelerate the search for the maximum exact match, POSIX-based multithreading is used for parallel processing. In addition, in order to shorten the execution time of the program, the FM index and the longest common prefix array are serialized. Experimental results show that the algorithm can greatly accelerate the search for the maximum exact match, and effectively accelerate the splicing alignment of RNA sequencing data.

Description

Technical Field

[0001] This invention relates to sequence alignment acceleration algorithms in bioinformatics, specifically a maximum exact match parallel search algorithm for accelerating the alignment of third-generation RNA sequencing data. Background Technology

[0002] With the emergence of third-generation sequencing technology, the field of genomics is undergoing a historic transformation. Unlike the first two generations of sequencing, third-generation sequencing generates longer reads that span the entire transcriptome, enabling researchers to decipher more complex gene structures and delve into intricate alternative splicing events. Simultaneously, the field of genomics is entering a post-genomic era, characterized by comprehensive genome sequencing to better understand the subtle differences in gene expression and regulation.

[0003] RNA sequencing has a profound impact on revealing the complexity of gene expression. The rise of the post-genomic era has provided a powerful impetus for uncovering complex gene regulatory mechanisms and has driven a profound restructuring of the dimensions of RNA sequencing data analysis. With in-depth research on alignment algorithms for large amounts of RNA sequencing data, current algorithms have made significant progress in accuracy, thereby promoting the widespread application of RNA sequencing and directly leading to a surge in sequencing data. This surge has ushered in the era of big data in genomics, posing new challenges to algorithm speed in order to process massive amounts of data faster and more efficiently. Modern computers typically possess multi-core capabilities; therefore, researching how to fully utilize CPUs to accelerate alignment algorithms is of significant practical importance for accelerating the analysis of RNA sequencing data. Summary of the Invention

[0004] The purpose of this invention is to accelerate the alignment of third-generation RNA sequencing data. To this end, a parallel maximum exact match search algorithm is provided. This algorithm features fast running speed, low memory consumption and high data throughput, which can effectively accelerate the alignment process and thus help the rapid analysis of RNA sequencing data.

[0005] The technical solution for realizing the present invention is as follows:

[0006] A maximum exact match parallel search algorithm for accelerating third-generation RNA sequencing data alignment includes the following steps:

[0007] (1) Collect datasets, including 3 real datasets and 1 simulated dataset. Each dataset includes a query sequence ref of the reference genome and third-generation RNA sequencing data reads;

[0008] (2) Perform BWT transformation on the reference genome sequence ref, where BWT is the Burrows-Wheeler Transform, which sorts all the rotation strings of the reference sequence to obtain an array of rotation strings, so that similar characters can be grouped together;

[0009] (3) Construct the FM index of the reference sequence, where the FM index refers to the Ferragina-Manzini index. The core is to construct the LF mapping based on the BWT rotating string array, which maps a certain position of the BWT rotating string array to the corresponding position of the original string.

[0010] (4) Calculate all suffixes of the reference sequence, sort the suffixes, and obtain the suffix array suffixes;

[0011] (5) Calculate the longest common prefix of adjacent suffixes in the suffix array to obtain the longest common prefix array LCP;

[0012] (6) Sample the longest common prefix array to obtain the sampled longest common prefix array SLCP;

[0013] (7) Serialize the longest common prefix array (SLCP) of the FM index and the samples to obtain the corresponding binary files index.bwt and slcpArray.bwt;

[0014] (8) Load a specified number of data from RNA sequencing data reads, wherein the number of data loaded in each batch is batchSize. When the remaining data volume is less than batchSize, all data are loaded.

[0015] (9) Initialize the multithreading threads[numThreads], set the parameters of each thread, including the start index and end index of the sequencing data to be processed by the thread, where numThreads is the number of threads used to find the maximum exact match;

[0016] (10) Create multiple threads, each thread finds the maximum exact match between the assigned sequencing data and the query sequence, where the maximum exact match refers to the longest consecutive equal substring between the RNA sequencing data and the reference sequence;

[0017] (11) Synchronize multiple threads and wait for all threads to finish executing;

[0018] (12) Merge the maximum exact match found by each thread into the output query result file;

[0019] (13) Repeat steps (8) to (12) until all RNA sequencing data in the dataset has been processed.

[0020] For the parallel search algorithm involved in this invention, runtime, speedup ratio, peak memory usage, and number of read segments processed per second are used as evaluation metrics.

[0021] The beneficial effects of this invention are: it provides a parallel maximum exact match search algorithm for accelerating the alignment of third-generation RNA sequencing data. This algorithm proposes a multi-threaded parallel maximum exact match search strategy, significantly reducing the time required to search for the maximum exact match between RNA sequencing data and the genome, while also optimizing the program's memory usage. Furthermore, by serializing the reference genome index used to find the maximum exact match, the time overhead of index creation is almost eliminated, thus enabling efficient handling of large-scale datasets. Applying this parallel maximum exact match search algorithm to the seed step of the RNA sequencing data splicing alignment algorithm significantly reduces the program's runtime without affecting the overall accuracy of the algorithm. Attached Figure Description

[0022] Figure 1 A flowchart of a maximum exact match parallel search algorithm for accelerating alignment of third-generation RNA sequencing data;

[0023] Figure 2 To improve the speedup of the parallel search algorithm for maximum exact match on different datasets with different batch sizes and different numbers of threads;

[0024] Figure 3 The number of segments per second processed by the splicing comparison algorithm on different datasets with different numbers of processes and threads. Detailed Implementation

[0025] To enable those skilled in the art to better understand the content of this invention, the invention will be further illustrated below through embodiments.

[0026] The following is an example:

[0027] A maximum exact match parallel search algorithm for accelerating third-generation RNA sequencing data alignment includes the following steps:

[0028] (1) Data sets were collected, including three real datasets (ALZ, DROS, and SIRV) and one simulated dataset (SIM). Each dataset included a reference genome query sequence (ref) and third-generation RNA sequencing reads. The ALZ dataset used Iso-Seq sequencing technology, the DROS and SIRV datasets used ONT sequencing technology, and the SIM dataset was obtained through simulated sequencing of human transcripts. The datasets were categorized based on their reference genome size: ALZ and SIM were large-scale datasets (3GB), DROS was a medium-scale dataset (140MB), and SIRV was a small-scale dataset (221KB).

[0029] (2) Perform a Burrows-Wheeler Transform (BWT) on the reference genome sequence `ref`. The BWT sorts all the rotation strings of the reference sequence to obtain an array of rotation strings, thus grouping similar characters together. For the query sequence `ref`, construct a matrix where each row represents the result of performing a cyclic rotation on `ref`, as shown in the following formula:

[0030] M(i,j)=ref[(i+j)mod len(ref)]

[0031] Where M(i,j) represents the element in the i-th row and j-th column of the matrix, and len(ref) represents the length of the query sequence ref.

[0032] Sort all rows of matrix M according to the lexicographical order of their last characters, as shown by the following formula:

[0033] M' = Sort(M[:,len(ref)-1])

[0034] Where M' represents the sorted matrix.

[0035] Extracting BWT can be represented by the following formula:

[0036] BWT = M'[:,len(ref)-1]

[0037] Where BWT is a string consisting of the last characters of the sorted matrix;

[0038] (3) Construct the FM index of the reference sequence, where the FM index refers to the Ferragina-Manzini index. The core is to construct the LF mapping based on the BWT rotating string array, which maps a certain position of the BWT rotating string array to the corresponding position of the original string.

[0039] (4) Calculate all suffixes of the reference sequence, sort the suffixes, and obtain the suffix array suffixes, where a suffix refers to any substring of ref that starts at a certain position, including the starting position and all characters thereafter;

[0040] (5) Calculate the longest common prefix of adjacent suffixes in the suffix array to obtain the longest common prefix array LCP, where LCP[i] represents the length of the longest common prefix between suffixes[i] and suffixes[i-1].

[0041] (6) Sample the longest common prefix array to obtain the sampled longest common prefix array SLCP;

[0042] (7) Serialize the longest common prefix array (SLCP) of the FM index and the samples to obtain the corresponding binary files index.bwt and slcpArray.bwt;

[0043] (8) Load a specified number of data from RNA sequencing data reads, wherein the number of data loaded in each batch is batchSize. When the remaining data volume is less than batchSize, all data are loaded.

[0044] (9) Initialize the multithreading threads[numThreads], set the parameters of each thread, including the start index and end index of the sequencing data to be processed by the thread, where numThreads is the number of threads used to find the maximum exact match;

[0045] (10) Create multiple threads, each thread finds the maximum exact match between the assigned sequencing data and the query sequence, where the maximum exact match refers to the longest consecutive equal substring between the RNA sequencing data and the reference sequence;

[0046] (11) Synchronize multiple threads and wait for all threads to finish executing;

[0047] (12) Merge the maximum exact match found by each thread into the output query result file;

[0048] Steps (8) to (12) can be represented by the following pseudocode:

[0049]

[0050] (13) Repeat steps (8) to (12) until all RNA sequencing data in the dataset has been processed;

[0051] (14) Evaluate the performance of the maximum exact match parallel search algorithm by running the algorithm on four datasets. For example, for the ALZ dataset, execute the following command:

[0052] . / PMEM-l 20refs_sequences.faALZ.fasta-oALZ-p 50000-t 16

[0053] Where "-l" is the minimum length for the maximum exact match, "-o" is the output query result file, "-p" is the batch size, and "-t" is the number of threads;

[0054] The peak runtime and memory usage of the search algorithm on various datasets were statistically analyzed, and the speedup and number of reads processed per second were calculated. For example, Figure 2 The figure shows the speedup of the parallel search algorithm for maximum exact match on different datasets with different batch sizes and different numbers of threads. As can be seen from the figure, the algorithm can effectively speed up the search for maximum exact match, and the speedup is more obvious on the two large-scale datasets, ALZ and SIM. That is, it has a significant advantage in speeding up the search for maximum exact match on large-scale datasets.

[0055] (15) Apply the maximum exact match parallel search algorithm to the seed step of the RNA sequencing data splicing alignment algorithm, and run the algorithm on four datasets. For example, for the ALZ dataset, execute the following command:

[0056] python uLTRA align Homo_sapiens.GRCh38.dna.primary_assembly.faALZ.fasta / ALZ-15--index / Homo_sapiens_GRCh38--prefix ALZ-15cores--isoseq--t15--min_mem 20--min_acc 0.8--dropoff 0.95--max_loc 5--disable_mm2--keep_temporary_files

[0057] Similarly, the runtime and peak memory usage of the splicing alignment algorithm were statistically analyzed on various datasets, and the speedup and number of reads processed per second were calculated. For example, Figure 3 The figure shows the number of reads per second processed by the splicing alignment algorithm on different datasets with different numbers of processes and threads. As can be seen from the figure, compared with the single-threaded maximum exact match parallel search algorithm can significantly improve the number of reads per second processed in the seed step of the splicing alignment algorithm, and the improvement is more obvious on the two large-scale datasets ALZ and SIM. That is, it has a significant advantage in accelerating splicing alignment on large-scale datasets.

Claims

1. A maximum exact match parallel search algorithm for accelerating the alignment of third-generation RNA sequencing data, characterized in that, Includes the following steps: (1) Collect datasets, including 3 real datasets and 1 simulated dataset. Each dataset includes a query sequence ref of the reference genome and third-generation RNA sequencing data reads; (2) Perform BWT transformation on the reference genome sequence ref. BWT is a Burrows-Wheeler Transform, which sorts all the rotation strings of the reference sequence to obtain an array of rotation strings, so that similar characters can be grouped together. (3) Construct the FM index of the reference sequence, where the FM index refers to the Ferragina-Manzini index. The core is to construct the LF mapping based on the BWT rotating string array, which maps a certain position of the BWT rotating string array to the corresponding position of the original string. (4) Calculate all suffixes of the reference sequence, sort the suffixes, and obtain the suffix array suffixes; (5) Calculate the longest common prefix of adjacent suffixes in the suffix array to obtain the longest common prefix array LCP; (6) Sample the longest common prefix array to obtain the sampled longest common prefix array SLCP; (7) Serialize the longest common prefix array (SLCP) of the FM index and the samples to obtain the corresponding binary files index.bwt and slcpArray.bwt; (8) Load a specified number of data from RNA sequencing data reads, wherein the number of data loaded in each batch is batchSize. When the remaining data volume is less than batchSize, all data are loaded. (9) Initialize the multithreading threads[numThreads], set the parameters of each thread, including the start index and end index of the sequencing data to be processed by the thread, where numThreads is the number of threads used to find the maximum exact match; (10) Create multiple threads, each thread finds the maximum exact match between the assigned sequencing data and the query sequence, where the maximum exact match refers to the longest consecutive equal substring between the RNA sequencing data and the reference sequence; (11) Synchronize multiple threads and wait for all threads to finish executing; (12) Merge the maximum exact match found by each thread into the output query result file; (13) Repeat steps (8) to (12) until all RNA sequencing data in the dataset has been processed.

Images