1. Field
The present description relates to aligning long sequences or patterns to find matches in sub-sequences or in portions and, in particular to using a grid cache and local start points to quickly find alignments of very long sequences.
2. Related Art
Sequence alignment is an important tool in signal processing, information technology, text processing, bioinformatics, acoustic signal and image matching, optimization problems, and data mining, among other applications. Sequence alignments may used to match sounds such as speech maps to reference maps, to match fingerprint patterns to those in a library and to match images against known objects. Sequence alignments may also be used to identify similar and divergent regions between DNA and protein sequences. From a biological point of view, matches point to gene sequences that perform similar functions, e.g. homology pairs and conserved regions, while mismatches may detect functional differences, e.g. SNP (Single Nucleotide Polymorphism).
Although efficient dynamic programming algorithms have been presented to solve this problem, the required space and time still pose a challenge for large scale sequence alignments. As computers become faster, longer sequences may be matched in less time. Multiple processor, multiple core, multiple threaded, and parallel array computing systems allow for still longer sequences to be matched. However, expanding uses of sequence alignment in information processing and other fields creates a demand for still more efficient algorithms. In bioinformatics, for example there is a great variety of organisms and millions of base pairs in each chromosome of most organisms.
The various advantages of the embodiments of the present invention will become apparent to one skilled in the art by reading the following specification and appended claims, and by referencing the following drawings, in which:
1. Introduction
In one embodiment, the invention for large-scale sequence alignment may be referred to as “SLSA” (Sequential Linear Space Algorithm). In SLSA, re-calculations are reduced by grid caches and global and local start points thereby improving overall performance. First, a whole similarity matrix H(i, j) is calculated in a linear space. The information on grids, including global and local start points and similarity values, are stored in grid caches. Then, the whole alignment problem is divided into several independent sub-problems. If a sub-problem is small enough, it will be solved directly. Otherwise, it will be further decomposed into several smaller sub-problems until the smaller sub-problems may be solved in the available memory. Using the global start points, several (k) near-optimal non-intersecting alignments between the two sequences can be found at the same time.
The grid cache and global and local start points used in SLSA, are efficient for large-scale sequence alignment. The local start points and grid cache divide the whole alignment problem into several smaller independent sub-problems, which dramatically reduces the re-computations in the backward phase and provides more potential parallelisms than other approaches. In addition, global start points allow many near-optimal alignments to be found at the same time without extra re-calculations.
In another embodiment, the invention for large-scale sequence alignment may be referred to as “Fast PLSA” (Fast Parallel Linear Space Alignment). Based on the grid cache and global and local start points, mentioned above, Fast PLSA provides a dynamic task decomposition and scheduling mechanism for parallel dynamic programming. Fast PLSA reduces sequential computing complexity by introducing the grid cache and global and local start points, and provides more parallelism and scalabilities with dynamic task decomposition and scheduling mechanisms.
Fast PLSA may be separated into two phases: a forward phase and a backward phase. The forward phase uses wave front parallelism to calculate the whole similarity matrix H(i, j) in linear space. The alignment problem may then be segmented into several independent sub-problems. The backward phase uses dynamic task decomposition and scheduling mechanisms to efficiently solve these sub-problems in parallel. This scheme can achieve automatic load balancing in the backward trace back period, tremendously improving the scalability performance especially for large scale sequence alignment problems.
2. Sequential LSA
Referring again to embodiments of the invention that may be characterized as Sequential LSA, for two sequences S1 and S2 with length l1 and l2, the Smith-Waterman algorithm (Temple F. Smith and Michael S. Waterman, Identification of Common Molecular Sub-sequences, Journal of Molecular Biology, 147:195-197 (1981)) computes a similarity matrix H(i, j) to identify optimal common sub-sequences by the using equation set 1, below:
where 1≦i≦l1, 1≦j≦l2 and sbt( ) is the substitution matrix of cost values. Affine gap costs are defined as follows: α is the cost of the first gap, and β is the cost of the following gaps. H(i, j) is the current optimal similarity value ending at position (i, j). E and F are the cost values from a vertical or horizontal gap respectively. An example is illustrated in
The memory required to solve the Smith-Waterman algorithm has been characterized as O(l1×l2), i.e. some independent factor (O) multiplied by the product of the length of each sequence (l1, l2). To align sequences with several hundred million elements, e.g. genome alignment, would lead to a memory requirement of several terabytes. Various different approaches have been developed to reduce the memory requirement. These usually increase the processing demands or reduce accuracy.
Fast LSA (Adrian Driga, Paul Lu, Jonathan Schaeffer, Duane Szafron, Kevin Charter and Ian Parsons, Fast LSA: A Fast, Linear-Space, Parallel and Sequential Algorithm for Sequence Alignment, In the International Conference on Parallel Processing, (2003)) uses some extra space called grid cache to save a few rows and columns of the similarity matrix H (see e.g.
2.1 k Near-Optimal Alignments and Global Start Points
The Smith-Waterman algorithm only computes the optimal local alignment result. However, the detection of near-optimal local alignments is particularly important and useful in practice. Global start point information may be used to find these different local alignments. The recurrences equation of the global start points is slightly inconvenient since it requires more computation and memory. However, the recurrences equation may be simplified as described below.
For each point (i, j) in the similarity matrix H, define the global start point Hst(i, j) as the starting point of the local alignment path ending at point (i, j). Similar to Eq (1), the values of Hst(i, j) may be calculated using the recurrence equations of equation set 2, below:
In order to determine k near-optimal alignments, the k highest similarity scores with different global start points are recorded during the forward period. If one of the k highest scores ends at a point (imax, jmax), then its global start point Hst(imax, jmax) can be easily obtained according to the stored information. The near optimal paths can be traced back in the rectangle of the two points. These near optimal paths with k highest scores will not intersect with each other. If the near optimal paths do intersect, then they must have the same global start point.
Using the start points, all of the k near-optimal alignments may be found at the same time without introducing extra re-computations. In addition, both the global and local alignment problem may be solved.
2.2 Grid Cache and Local Start Points.
For many processing systems, the system memory is not large enough to contain the complete similarity matrix for long sequences. A partial similarity matrix H may then be re-computed in the backward phase. To reduce re-calculations, a few columns and rows of the matrix H may be stored in k×k grid caches.
The sub-problems can be processed only after the last sub-problem, the adjacent bottom-right grid cache 222, is solved. After all of the sub-problems are solved recursively, the sub-paths may be concatenated to form the full optimal alignment path 222.
In combination with grid caches, local start points may be used to generate smaller and independent sub-problems. Similar to the global start point described above, the local start point of one point (i, j) may be defined as the starting position in its left/up grid of the local alignment ending at point (i, j). The local start point may be calculated by Eq (2) with different initialization on the grids. Using the grid cache and local start points, the whole alignment problem can be divided into several independent sub-problems.
As shown in
2.3 Solving Sub-Problems
In order to improve the trade-off between time and space, a block may be used as the basic matrix filling and tracing path unit. The block, similar to a 2D matrix, denotes a memory buffer which is available for solving small sequence alignment problems. If a problem or sub-problem is small enough, it may be directly solved within a block. Otherwise it will be further decomposed into several smaller sub-problems until the sub-problems are small enough to easily be solved. Since the start and end points are fixed in the sub-problem, it becomes a global alignment problem. For global alignment, the computation of the score H(i, j) may be given by the recurrences of equation set 3, below.
In order to improve performance, the block size may be tuned to suit different memory size and cache size configurations. All of the sub-problems may be solved in parallel for faster speed since they are independent of each other. After all the sub-problems are solved, the traced sub-paths may be concatenated to produce a full optimal alignment path.
In sum, Sequential LSA, as described herein represents a fast linear algorithm for large scale sequence alignment. The joint contribution of the grid cache and global and local start points, allow a large-scale alignment problem to be recursively divided into several independent sub-problems until each independent sub-problem is small enough to be solved. This approach dramatically reduces the re-computations in the backward phase and provides more parallelism. In addition, using global start points can efficiently find k near-optimal alignments at the same time.
2.4 Pseudo-Code for Sequential LSA
The Sequential LSA approach described above may be represented, in one example, by the flow chart of
The forward phase begins with block 312 by calculating a similarity matrix H for an input pair of sequences that are to be compared. Information about the grids of matrix H is then stored in grid caches at block 314. This information may include global and local start points and similarity values. The grids may be based on a block size and grid division that may be set before the alignment process is started. At block 316, the ending point that has the maximum score in H is found and then the optimal path may be identified based on the global starting point and the found ending point. The local points may also be found based on the optimal path at block 318.
The whole problem may then be divided into sub-problems using the local start points at block 320. The sub-problems may then be pushed into problem queues at block 322. If the sub-problem can be solved within a block size, then it is pushed to a solvable problem queue. If the sub-problem cannot be solved within a block size, then it is pushed into an unsolvable problem queue.
The backward phase begins by processing the problems in the unsolvable problem queues at block 324. This may be done in the same way as in the forward process. The processing may divide the unsolvable problems into smaller sub-problems based on local starting points until they can be solved within a block size. Then the problems may be pushed into solvable problem queues. At block 326, the sub-problems in the solvable problem queues are solved and the sub-paths are traced backwards to find the sub-alignment paths. At block 328, when all the sub-problems are solved then all the sub-alignments are concatenated into a final alignment solution.
A process such as that of
Input: sequence 1, 2; block size (h, w); grid division (rows, cols); Output: optimal alignment path
Initialize unsolvable problem queue and solvable problem queue to empty.
1. Forward process:
1.1 Calculate the whole similarity matrix H in linear space. The information on grids, including global/local start points and similarity value, are stored in the grid caches.
1.2 Find the ending point with max score in H and get the optimal path's global/local points from the ending point.
1.3 Divide the whole problem into independent sub-problems by these local start points
1.4 Push these sub-problems into a queue depending on whether they can be directly solved within a block size or not.
2. Backward process:
2.1 Process each unsolvable sub-problem in the unsolvable problem queue using the same strategy as the forward process until the sub-problems are solvable.
2.2 Solve the sub-problems in the solvable problem queue to trace back the sub-paths.
2.3 If all the sub-problems are solved, concatenate all solutions into output alignment path.
3. Fast Parallel Linear Space Algorithm (Fast PLSA)
In another embodiment, an approach denoted FastPLSA uses the grid cache and global start point described above to reduce the sequential execution time and to provide more parallelism especially when coping with the trace back phase. In addition, Fast PLSA can output several near-optimal alignment paths after one full matrix filling process. It also introduces several tunable parameters so that the whole process can adapt to different hardware configurations, such as cache size, main memory size and the different communication interconnections, data rates and speeds. FastPLSA is able to use more available memory to reduce the re-computation time for the trace back phase.
To further improve alignment performance, the Sequential LSA algorithm may be parallelized using a Fast PLSA approach. Large scale sequence alignments may be mapped to a parallel processing architecture in two parts. The first part is forward calculating the whole similarity matrix and the second part is backward solving sub-problems to find trace path phase.
3.1 Forward Phase
In the forward phase, block 410 of
The forward phase begins with initializing all of the values for memory grid size, problem queues etc. at block 412. The computation of each block follows the dynamic programming of equation set 1 to fill the block matrix. The whole similarity matrix may be built by first initializing the values of the leftmost column and topmost row at block 414. The top left block may then be computed immediately. Since each block has dependencies on its adjacent left, upper left and upper blocks according to equation set 1, it may be processed at block 416 after its adjacent related blocks are computed. Based on such a dependency model, a wave front communication pattern can be used in the parallelization of the similarity matrix.
The wave front moves in anti-diagonals as depicted in FIGS. 5A, 5B,and 5C.
The wave front computation may be parallelized in several different ways depending upon the particular parallel processing architecture that will be used. On fine-grained architectures such as shared memory systems, the computation of each cell or a relatively smaller block within an anti-diagonal may be parallelized. This approach works better for very fast inter-processor communications since the granularity for each processing unit is extremely small. On the other hand, for distributed memory systems such as PC clusters, it may be more efficient to assign a relatively larger block to each processor. In one example, two parameters h and w are used to denote the height and width of each block in terms of cells. These may be tuned to adapt to different architectures.
In the Fast PLSA example of
Referring to
The second processor may use the transferred margin as the initial top margin of block (1,0) and the right margin of block (0,0) may be used as the initial left margin of block (0,1). In this way, block (0,1) and block (1,0) may be computed by the first and second processors simultaneously as soon as block (0,0) is completed.
Similarly, additional processors may process additional blocks at the same time. The processing of these tiles advances on a diagonal wave front 516, 517, 518, and more of the processors can be added to work in parallel as the diagonal wave front progresses. If there are P processors in total and it requires one time step per tile, then all P processors are operating after (P-1) time steps.
Along with the block computing, the grid cache may be saved when a part of the grid columns or rows is within a computing block. Since the grid cache is distributed among all the processors, a procedure denoted in
In many implementations, each block will have two communication operations, receiving the bottom marginal data from the upper block and sending the upper block's marginal data to the bottom block. The communication overhead may be reduced especially in a PC cluster by using non-blocking receive message passing operations to overlap the communication overhead with computing. The receive message passing operations may work like a pipeline block by block until the whole similarity matrix H is filled. This minimizes the communication cost and delivers better parallelization performance.
3.2 Backward Phase
After the forward phase 410, a series of independent sub-problems are stored in the unsolvable problem queue. For each sub-problem, a global alignment technique may be used in a backward phase 430 to solve these sub-problems and concatenate these alignments to the optimal path. Sub-problem alignment may be solved by a repeated process of the forward phase. However, there may be several differences between forward and backward phases.
a) Since the start and end points of the optimal alignment paths are unknown in the forward phase, a Smith-Waterman algorithm may be used to fill the whole similarity matrix and find all the sub-problems. In the backward phase, each sub-problem has fixed start and end points, so, for example, a Needleman-Wunsch algorithm may be used to find global alignment of these sub-problems. Saul B. Needleman and Christian D. Wunsch. “A General Method Applicable to the Search for Similarities in the amino acid Sequence of Two Sequences” Journal of Molecular Biology, 48:443-453 (1970).
b) Different parallel schemes may be used in the forward phase and the backward phase. The forward phase may use wave front parallelism as described above. In the backward phase, since all the sub-problems are independent of each other, more factors, may be considered such as the size of sub-problems and the number of processors. Factors may also be combined to derive better parallel schemes. Attention to the load balance performance to efficiently use all the processors in the backward period may be particularly effective because the backward phase's granularity is much finer than the forward phase's granularity.
c) For large scale alignment problems, in general, the problem may be divided into several sub-problems in the forward phase. In the backward phase, if the sub-problem size is smaller than the block size, it may be directly solved by using the full matrix filling method. Otherwise, approaches similar to those used in the forward phase may be used to subdivide sub-problems.
The differences between the forward phase and the backward phase allow the two phases to be tailored differently to improve computational efficiency, accuracy and speed. In one implementation, the sub-problems may be evenly and independently distributed to all of the processors. Each processor then works on a sub-problem using the sequential methods described above. After the sub-problems are solved, the processors collect the sub-alignments together and concatenate them to the optimal alignment.
To better balance the processing load among the processors, each sub-problem may first be recursively decomposed in a wave front parallel scheme until all the descendant sub-problems are reduced to the block size and can be quickly solved. This recursive decomposition may be applied to each sub-problem in turn. This scheme is particularly effective for small scale processors or large scale sub-problems. Many modifications and variations may be made to these and the other approaches described above to consider both the load balance as well as the granularity of the problems in the backward parallel phase, and to design a flexible scheme to partition tasks equally for all the processors.
In one embodiment as shown in
The “balanced state” means that all of the sub-problems may be distributed roughly equally to all the processors within some threshold (e.g. 20%). In other words, the “balanced state” indicates that the difference of the sum area of the sub-problems assigned to each processor are within the threshold value. If, for example, the unsolved sub-problem queue consists of four sub-problems of different sizes (100×100, 50×50, 70×70 and 110×100) to be assigned to two different processors, then to evenly distribute tasks between these two processors, the first processor may be assigned the 100×100, and 70×70 tasks, and the second processor may be assigned the 50×50 and the 110×100 sub-problems respectively. The size difference ratio may be computed for the two processors, and the value (14900-13500)/13500=10.3% is smaller than the default threshold. Therefore, the unsolved sub-problem queue is in the “balanced state”.
In one embodiment, a formula may be applied to determine whether the sub-problems in the queue are in the “balanced state” as shown in equation set 4, below:
Sizeaverage=ΣSizei/M, where the sum is for i=0 to M (4)
|(Sizepj−Sizeaverage)/Sizeaverage|<Threshold, 1≦j≦M
where M is the total processor number and N is the sub-problem number. Sizei is the area for each sub-problem and Sizepj is the total area of sub-problems assigned to the jth processor. If the difference of the sum area of the sub-problems assigned to each processor is within the Threshold value (in one embodiment, a default value may be 20%), the sub-problems can be considered to enter the “balanced state”, indicating that the sub-problems are distributed equally to each processor.
If the unsolved sub-problem queue is not in the “balanced state”, then the largest size sub-problem from the queue may be found and decomposed into several smaller descendant sub-problems with wave front parallelism. After that, the descendant sub-problems may be pushed back into the unsolved problem queue. The balanced state test may then be iterated to detect whether the queue is again in the “balanced state” or not.
Referring to
In
After the unsolved sub-problem queue is in the “balanced state,” the individual solving sub-problem phase 434 of
3.4 Pseudocode for Fast PLSA
A process such as that of
Input: sequence 1, 2; block size (h, w); grid division (rows, cols); Output: optimal alignment path
Forward process:
1.1 Calculate the whole similarity matrix H in linear space with wave front parallel scheme.
1.2 The information on grids, including global/local start points and similarity value, are stored in the grid caches.
1.3 Collect all the distributed grid cache information to the root processor.
1.4 Find the ending point with max score in H and get the optimal path's global/local start points from the ending point.
1.5 Divide the whole problem into independent sub-problems by these local start points
1.6 Push all these sub-problems into the “unsolved queue”
2. Backward Process:
2.1 If all the sub-problems in the “unsolved queue” can be distributed to the processors equally, pick out the largest sub-problem and subdivide it into a series of smaller sub-problems using the same strategy as the forward process.
2.2 Push all of those decomposed sub-problems back into the “unsolved queue”, go back to 2.1
2.3 Otherwise, go directly into the individual work phase, where all the sub-problems in this queue will be proximately assigned to the working processors.
2.4 Each processor will work independently to find the sub alignment paths for the assigned sub-problems.
3. Concatenate all the sub alignments individually on each processor, and finally, merge them together into the final one.
The Fast PLSA approach produces k near-optimal maximal non-intersecting alignments within one forward and one backward phase. The speedup in k alignments (k>1) is usually better than for a single alignment. This may be because the forward phase execution time is relatively stable and more sub-problems can be generated when the number of output alignments is increased. In the example of
The described approaches allow for long sequence alignments to be performed more quickly using linear space. A trade is made to increase space to reduce time. The local start points and grid cache can divide the whole sequence alignment problem into several independent sub-problems, which dramatically reduces the re-computations of the trace back phase and provides more parallelism. The dynamic task decomposition and scheduling mechanism can efficiently solve the sub-problems in the backward phase. This tremendously improves the scalability performance and minimizes the load imbalance problem especially for large scale sequence alignment.
4 Processing Environment
The approaches described above may be carried out on a variety of different processing environments. In one embodiment, a 16-node PC cluster interconnected with a 100 Mbps Ethernet switch may be used. Each node has a 3.0 GHz Intel Pentium-4 processor with 512 KB second-level cache and 1 GB memory. The RedHat 9.0 Linux operation system and MPICH-1.2.5 message passing library (Message Passing Interface from Mathematics and Computer Science Division, Argonne National Laboratory, Illinois) may be used as the software environment. The sequence alignment routines may be written in C++ or any other programming language or implemented in specialized hardware.
The particular architecture of
The MCH may also have an interface, such as a PCI Express, or AGP (accelerated graphics port) interface to couple with a graphics controller 341 which, in turn, provides graphics and possible audio to a display 337. The PCI Express interface may also be used to couple to other high speed devices. In the example of
The ICH 365 offers possible connectivity to a wide range of different devices. Well-established conventions and protocols may be used for these connections. The connections may include a LAN (Local Area Network) port 369, a USB hub 371, and a local BIOS (Basic Input/Output System) flash memory 373. A SIO (Super Input/Output) port 375 may provide connectivity to a keyboard, a mouse, and other I/O devices. The ICH may also provide an IDE (Integrated Device Electronics) bus or SATA (serial advanced technology attachment) bus for connections to disk drives 387, or other large memory devices.
The particular nature of any attached devices may be adapted to the intended use of the device. Any one or more of the devices, buses, or interconnects may be eliminated from this system and others may be added. For example, video may be provided on a PCI bus, on an AGP bus, through the PCI Express bus or through an integrated graphics portion of the host controller.
5. General Matters
A lesser or more equipped optimization, process flow, or computer system than the examples described above may be preferred for certain implementations. Therefore, the configuration and ordering of the examples provided above may vary from implementation to implementation depending upon numerous factors, such as the hardware application, price constraints, performance requirements, technological improvements, or other circumstances. Embodiments of the present invention may also be adapted to other types of data flow and software languages than the examples described herein. The methods described above may be implemented using discrete hardware components or as software.
Embodiments of the present invention may be provided as a computer program product which may include a machine-readable medium having stored thereon instructions which may be used to program a general purpose computer, mode distribution logic, memory controller or other electronic devices to perform a process. The machine-readable medium may include, but is not limited to, floppy diskettes, optical disks, CD-ROMs, and magneto-optical disks, ROMs, RAMs, EPROMs, EEPROMs, magnet or optical cards, flash memory, or other types of media or machine-readable medium suitable for storing electronic instructions. Moreover, embodiments of the present invention may also be downloaded as a computer program product, wherein the program may be transferred from a remote computer or controller to a requesting computer or controller by way of data signals embodied in a carrier wave or other propagation medium via a communication link (e.g., a modem or network connection).
In the description above, numerous specific details are set forth. However, it is understood that embodiments of the invention may be practiced without these specific details. For example, well-known equivalent components and elements may be substituted in place of those described herein, and similarly, well-known equivalent techniques may be substituted in place of the particular techniques disclosed. In other instances, well-known circuits, structures and techniques have not been shown in detail to avoid obscuring the understanding of this description.
While the embodiments of the invention have been described in terms of several examples, those skilled in the art may recognize that the invention is not limited to the embodiments described, but may be practiced with modification and alteration within the spirit and scope of the appended claims. The description is thus to be regarded as illustrative instead of limiting.