JP2012032975A - Optimal alignment computation device and program - Google Patents

Abstract

PROBLEM TO BE SOLVED: To accelerate dynamic programming (DP) alignment processing by parallelization with one processor, utilizing bit parallelism of word operations in a computer, without imposing restrictions on the lead length.SOLUTION: A range to compute a DP matrix is restricted to a near-diagonal region with the width corresponding to the word length in a computer, and information inside this region is divided into and represented by bits aligned in the direction orthogonal to the diagonal to achieve bit parallelization of computation for information in these bits.

Description

  The present invention speeds up the process of comparing and collating character string data such as character strings and DNA base sequences (including arbitrary data processed into a form that can be handled as character string data) with search character strings. Related to technology.

  In recent years, a massively parallel DNA sequencer has emerged that can read tens of millions of sequences in a single run based on a principle different from that of conventional capillary DNA sequencers. Was drastically reduced (Non-patent Document 1). Therefore, there is an increasing demand for a technique for analyzing a large amount of array data at high speed. Among them, genome mapping analysis is one of the most basic sequence data analysis. The base sequence (read) newly read by the sequencer is compared with the reference genome sequence in the database, and the position in the genome is compared. And mutations (base substitution / insertion / deletion) contained therein are identified.

  Since the read length of the initial massively parallel DNA sequencer was relatively short, about 30 bases, high-speed and high-accuracy genome mapping analysis was performed using a hash method or Burrows-Wheeler transformation (Non-patent Document 2).

  However, when the read length increases to 100 bases or more, a method such as the seed-and-extend strategy that has been used in genome mapping analysis for conventional DNA sequencers becomes necessary (non-patent literature). 3). In this method, in the first step, the position where the short sequence (seed) contained in the read appears in the genome is searched at high speed, and in the second step, the entire read is similar to the genome sequence around the seed appearance position. (That is, whether the correspondence between the seed and the genome can be extended to the entire read while allowing some variation). In this first step, since the seed length is short, a high-speed search method using a hash method, Burrows-Wheeler transform, etc. can be used. You can narrow down. In the second step, the genome partial sequence (target sequence) at the candidate position narrowed down in the first step is aligned with the read (query sequence) using dynamic programming (DP method), etc. Even when the lead length is long, the base mutation can be examined with high accuracy (Non-patent Document 8).

  However, in general, dynamic programming requires a large calculation cost proportional to the square of the array length. Therefore, speeding up the alignment process by dynamic programming is beneficial for speeding up the entire genome mapping analysis process.

  In order to speed up alignment processing by dynamic programming, a method of parallelizing processing using a plurality of processors is known (Non-Patent Document 4). At this time, it is necessary to consider the dependency relationship of data passed between the processors and to adopt a load balancing method that suppresses the inter-processor communication cost. Therefore, if a large number of processors are used, high parallelism can be obtained.

  By the way, in the current mainstream digital computer, an operation on an 8-byte integer (word) is a basic operation, but this can also be regarded as an 8-byte or 64-bit parallel operation. G. Myers showed how to use the parallelism of integer programming to parallelize the alignment process by dynamic programming with one processor and calculate in time proportional to the array length (Non-patent Document 5). .

  In G. Myers' method, in order to obtain an alignment between a query sequence (read) of length m and a target sequence (genome partial sequence) of length n, information on a DP matrix of size m × n is used as a column of length m. Processes related to the comparison of one character of the newly read target sequence and all the m characters of the query sequence are performed in m bits in parallel and expressed in the row direction (target sequence of the target sequence). Repeat this process while moving one base at a time).

This will be described in detail using equations. All the substitution, insertion, and deletion costs of one base are 1, and the DP matrix's (i, j) component C [i, j] (query sequence q = q1 q2 ... qm length i partial character For the optimal alignment cost of the sequence q1 q2 ... qi and the target sequence t = t1 t2 ... tn with substring t1 t2 ... tj of length j)
A Boolean variable indicating whether or not the equation C [i, j] −C [i, j−1] = 1 is satisfied, Ph [i, j],
A Boolean variable indicating whether or not the equation C [i, j] −C [i, j−1] = − 1 holds is represented by Mh [i, j],
A Boolean variable indicating whether or not the equation C [i, j] −C [i−1, j] = 1 is satisfied is Pv [i, j],
A Boolean variable indicating whether or not the equation C [i, j] −C [i−1, j] = − 1 holds is expressed as Mv [i, j],
Eq (i, j) is a Boolean variable that indicates whether the character (i-th character) qi of the index i of the query sequence matches the character tj of the index j of the target sequence
And a word variable that can be arranged in ascending order from the least significant bit for i = 1, 2, ..., m for these Boolean variables represented by 1 bit (false or true, ie 0 or 1),
(Equation 1)
Ph (j) = [Ph [m, j], ..., Ph [2, j], Ph [1, j]]
Mh (j) = [Mh [m, j], ..., Mh [2, j], Mh [1, j]]
Pv (j) = [Pv [m, j], ..., Pv [2, j], Pv [1, j]]
Mv (j) = [Mv [m, j], ..., Mv [2, j], Mv [1, j]]
Eq (j) = [Eq [m, j], ..., Eq [2, j], Eq [1, j]]
In order to simplify the symbol, j is fixed to one value, and the values of Ph (j), Mh (j), Pv (j), Mv (j), Eq (j) are respectively Abbreviated as Ph, Mh, Pv, Mv, Eq, Xv, Xh as auxiliary word variables,
(Equation 2)
Xh = ((Eq & Pv) + Pv) ^ Pv | Eq
Ph = Mv | ~ (Xh | Pv)
Mh = Pv & Xh
Ph << = 1
Mh << = 1
Xv = Eq | Mv
Pv = Mh | ~ (Xv | Ph)
Mv = Ph & Xv
Executes bitwise integer arithmetic and overwrites the values of Ph, Mh, Pv, and Mv. Here, for word variables a and b, a | b is a bitwise OR operation, a & b is a bitwise AND operation, a ^ b is a bitwise exclusive OR operation, and a + b is an integer Addition, ~ a represents a bitwise NOT operation, and a << = 1 represents a 1-bit left shift operation. Thereby, the values of Ph, Mh, Pv, and Mv are updated to the values of Ph (j + 1), Mh (j + 1), Pv (j + 1), and Mv (j + 1). This is the core part of aligning the alignment process bitwise, which speeds up the process.

  However, this method can be used directly only when the read length (query sequence length) m is less than or equal to the computer word length (64). In order to extend this method when the read length exceeds the word length of the computer, it is necessary to divide the read (query sequence) into a plurality of fragment sequences having a length equal to or less than the word length (Non-Patent Document 5, Non-Patent Document 5, Non-patent Document 5). Patent Document 6). At this time, the calculations for these fragment sequences are dependent on each other, and restrictions are imposed upon parallelization, resulting in a reduction in processing efficiency.

Rajendrani Mukhopadhyay.DNA sequencers: the next generation.Analytical Chemistry, Vol. 81, No. 5, March 1, 2009 Kouichi Kimura, Yutaka Suzuki, Sumio Sugano, Asako Koike.Computation of Rank and Select Functions on Hierarchical Binary String and Its Application to Genome Mapping Problems for Short-Read DNA Sequences.Journal of Computational Biology.November 2009, 16 (11): 1601 -1613. Wu TD, Watanabe CK. GMAP: a genomic mapping and alignment program for mRNA and EST sequences. Bioinformatics. 2005 May 1; 21 (9): 1859-75. Martins WS, Del Cuvillo JB, Useche FJ, Theobald KB, Gao GR.A multithreaded parallel implementation of a dynamic programming algorithm for sequence comparison.Pac Symp Biocomput. 2001: 311-22. G. Myers.A fast bit-vector algorithm for approximate string matching based on dynamic programming.Journal of the ACM, 46 (3): 395-415, 1999. HEIKKI HYYRO. BIT-PARALLEL COMPUTATION OF LOCAL SIMILARITY SCORE MATRICES WITH UNITARY WEIGHTS. International Journal of Foundations of Computer Science, 17 (6): 1325-1344, 2006. Kevin Judd McKernan, Heather E. Peckham, Gina L. Costa, et al. Sequence and structural variation in a human genome uncovered by short-read, massively parallel ligation sequencing using two-base encoding.Genome Res. 2009 19: 1527-1541 . Richard Durbin, Sean R. Eddy, Andrew Krogh, Graeme Mitchison, Biological Seuqnece Analysis, probabilistic models of proteins and nucleic acids, Cambridge University Press 1988 Sequence analysis, medical publication, 2001)

  The problem to be solved by the present invention is to make use of the bit parallelism of the word operation of a computer without setting a restriction on the read length and without using a method of reducing the processing efficiency such as dividing the read array. Another object is to provide a method for parallelizing and speeding up alignment processing by dynamic programming with one processor.

  In the method according to the present invention, in order to obtain an alignment between a query sequence (read) of length m and a target sequence (genome partial sequence) of length n, the computation of a DP matrix of size m × n is performed diagonally (row index i And the column index j are equal to each other), and the information in that area is expressed by disassembling the information in the direction perpendicular to the diagonal (opposite diagonal direction), and the bit of the word operation of the computer Using parallelism, the processing related to character comparison between the query sequence and the target sequence on the opposite corner is parallelized by one processor, and this processing is repeated while moving in the diagonal direction.

This will be described in detail using equations. C [i, j], Ph [i, j], Mh [i, j], Pv [i, j], Mv [i, j], Eq [i, j] Word variables that can be arranged in the opposite corner direction around
(Equation 3)
Ph (k) = [Ph [i, j] | i + j = k, −m <i−j ≦ m]
Mh (k) = [Mh [i, j] | i + j = k, −m <i−j ≦ m]
Pv (k) = [Pv [i, j] | i + j = k, −m <i−j ≦ m]
Mv (k) = [Mv [i, j] | i + j = k, −m <i−j ≦ m]
Eq (k) = [Eq [i, j] | i + j = k, −m <i−j ≦ m]
And Here, the order in which the bit variables are packed into the word variables is ascending order from the lower bits for j. Furthermore, to simplify the symbol, k is fixed to one value, and the values of Ph (k), Mh (k), Pv (k), Mv (k), Eq (k) are set to Ph, Abbreviated as Mh, Pv, Mv, Eq, Xv, Xh as auxiliary word variables,
When k is an even number,
(Equation 4)
Xv = Eq | Mv; Pv '= Mh | ~ (Xv | Ph); Mv' = Ph &Xv;
Xh = Eq | Mh; Ph '= Mv | ~ (Xh | Pv); Mh' = Pv &Xh;
When k is odd,
(Equation 5)
Ph >> = 1; Mh >> = 1;
Xv = Eq | Mv; Pv '= Mh | ~ (Xv | Ph); Mv' = Ph &Xv;
Xh = Eq | Mh; Ph '= Mv | ~ (Xh | Pv); Mh' = Pv &Xh;
Pv '<< = 1; Mv'<< = 1;
According to the above, a bit-parallel integer operation is performed, and the values of Ph ', Mh', Pv, and 'Mv' are calculated. These give the values Ph (k + 1), Mh (k + 1), Pv (k + 1), Mv (k + 1). This is the core part that makes the alignment process bit parallel.

  Similar to G. Myers' method, the alignment process is bit-parallelized and accelerated by one processor. The difference from G. Myers' method is that the bit parallel computation of dynamic programming is limited to the area around the diagonal of the DP matrix and moved in the diagonal direction, so there is no restriction on the read length, so the read length is equal to the word When the length is exceeded, there is no need to divide the lead. In addition, dynamic programming calculations for sequences with high similarity need only be performed around the diagonal of the DP matrix, so unnecessary calculations in areas away from the diagonal can be avoided, and processing efficiency is improved. improves. In addition, since this is a new bit parallelization method that can perform calculations in a symmetric manner with respect to rows and columns, in (Equation 4) and (Equation 5), the equation for the auxiliary variable Xh is symmetric with the equation for Xv. This can be calculated with a smaller number of machine instructions than the slightly complicated calculation formula of the auxiliary variable Xh in (Equation 2) by G. Myers.

In the optimal alignment calculation method by this invention, it is explanatory drawing which showed the application location of bit parallel calculation (Example 1). (Example 1) which is a figure which shows the structural example of the system which performs the high-speed comparison collation with the sequence data of a next-generation type DNA sequencer, and a reference genome sequence database. It is explanatory drawing which shows the positional relationship of the Boolean variable of the input side and output side in the cell periphery in DP matrix. It is explanatory drawing which illustrates the correspondence of the double index [i, j] whose sum is k, and the index of a word variable. It is explanatory drawing which illustrates the correspondence of the index of the Boolean variable Ph, Mh, Pv, Mv and the index of a word variable with respect to fixed k. It is explanatory drawing which shows that a difference arises according to whether the index of the Boolean variables Ph, Mh, Pv, and Mv around a cell is k even or odd. It is explanatory drawing which shows the method of calculating the Boolean variable Eq in bit parallel. It is explanatory drawing which shows the Boolean variable value preset in the exterior of a diagonal periphery area | region. It is explanatory drawing which shows the accumulation method of the difference value of DP cost for calculating | requiring the absolute value of DP cost.

[Example 1]
Hereinafter, embodiments of the present invention will be described in detail with reference to the drawings.
In this example, when a read (query sequence) and reference genome sequence data are given, a position (mapping position) in the genome sequence where a partial sequence similar to the read exists is obtained, and the genome around the query sequence and the mapping position is obtained. A method for obtaining the optimum alignment with the partial sequence (target sequence) will be described. Usually, since the lead sequence obtained from a DNA sequencer is a base sequence, it is assumed hereinafter that the lead sequence is a base sequence. As a special case, in a DNA sequencer adopting the 2-base encoding method, the lead sequence is a color sequence in which two adjacent bases are encoded in four colors (Non-patent Document 7). In such a case, it is possible to convert the reference genome sequence into a color sequence and apply a method similar to the method described in this embodiment by simple analogy. In addition, a method similar to the method described in this embodiment can be applied to the problem of searching a reference protein sequence database using an amino acid sequence as a query sequence, by simple analogy. In addition, a method similar to the method described in the present embodiment can be applied to a task of searching a document or other database including a character string or an image by simple analogy.

  FIG. 2 is a diagram illustrating a system configuration example of the present embodiment. Reference numeral 201 denotes a computer in which this system is operated, reads information of the reference genome sequence data 203 stored in the external storage device 202 into the main memory 209, and analyzes the DNA sample 205 using the DNA sequencer 204. The base sequence data (query sequence, read sequence) 206 obtained in this way is read into the main memory 209 via the input device 207, and processing parameters specified by the user (upper limit value M of alignment cost of similar sequences) ) Is read into the main memory 209 via the input device 208.

  The computer 201 performs a seed search process 210 on the read sequence and the reference genome sequence stored in the main memory 209 to obtain a position where a short sequence (seed sequence) included in the read appears in the genome sequence. As the seed sequence, it is sufficient to use a short partial sequence of about 12 to 20 bases in length from the lead sequence, and Burrows-Wheeler transformation or the like was used to determine the appearance position of the seed sequence. A known high-speed search method may be used (Non-Patent Document 2).

  The detailed comparison processing unit 211 performs an optimal alignment calculation 212 around the position where the seed sequence appears in the genome sequence. If the cost value of the alignment is equal to or less than the upper limit value M specified by the user, the lead sequence is included there. It is determined that there is a sequence similar to, and the output position information and alignment information in the genome sequence are transmitted to the output processing unit 213. The output processing unit 213 outputs the appearance position information (genome mapping position information) 215 in the genome sequence and the optimum alignment result 216 as search results to the output device 214.

  The cost of alignment is calculated as the sum of the costs of substitution, insertion, and deletion of one base included in the alignment as 1. If a seed sequence of length s starting from the position of the p base of the query sequence appears in the coordinate interval [x, x + s−1] in the genome sequence, the appearance position of the lead sequence is the coordinate interval [ x−p + 1−M, x−p + m + M]. Here, m represents the query sequence length, and M represents the upper limit of the alignment cost. The reason why the range of the appearance position of the lead sequence is expanded by ± M is that there is a possibility that insertion / deletion of a maximum of M bases is included before and after the appearance position of the seed sequence. In the optimal alignment calculation 212, a semi-global alignment is performed between the query sequence and the target sequence using the genome partial sequence in the coordinate interval [x−p + 1−M, x−p + m + M] as the target sequence. I do. That is, a correspondence relationship is determined between the entire query sequence and a partial sequence of the target sequence while allowing insertion / deletion. Specifically, the calculation is performed by calculating a DP matrix (Non-patent Document 8).

The query sequence length is m, and the target sequence length is n. For 1 ≤ i ≤ m, 1 ≤ j ≤ n, query string q = q1 q2 ... qm of length i partial string q1 q2 ... qi and target array t = t1 t2 ... tn According to G. Myers, the first (i, j) component of the DP matrix representing the optimal alignment cost with the substring t1 t2 ... tj of length j is C [i, j] (non-patent literature) 5),
A Boolean variable indicating whether or not the equation C [i, j] −C [i, j−1] = 1 is satisfied, Ph [i, j],
A Boolean variable indicating whether or not the equation C [i, j] −C [i, j−1] = − 1 holds is represented by Mh [i, j],
A Boolean variable indicating whether or not the equation C [i, j] −C [i−1, j] = 1 is satisfied is Pv [i, j],
A Boolean variable indicating whether or not the equation C [i, j] −C [i−1, j] = − 1 holds is expressed as Mv [i, j],
Eq (i, j) is a Boolean variable that indicates whether or not the character (i-th character) qi of the index i of the query sequence matches the character j of the index j (j-th character) tj of the target sequence.
It expresses. Further, C [i, 0] = i, C [0, j] = 0 are defined as boundary conditions corresponding to the semi-global alignment. For any 1 ≤ i ≤ m, 1 ≤ j ≤ n,
(Equation 6)
Mh = Mh [i−1, j]; Ph = Ph [i−1, j];
Mv = Mh [i, j−1]; Pv = Pv [i, j−1];
Mh '= Mh [i, j]; Ph' = Ph [i, j];
Mv '= Mv [i, j]; Pv' = Pv [i, j];
Eq = Eq [i, j];
For short, between these Boolean variables,
(Equation 7)
Xv = Eq or Mv; Pv '= Mh or not (Xv or Ph); Mv' = Ph and Xv;
Xh = Eq or Mh; Ph '= Mv or not (Xh or Pv); Mh' = Pv and Xh;
Is established (Non-patent Document 5). Here, Xv and Xh represent auxiliary Boolean variables. By using this relational expression repeatedly, all values of C [i, j] can be calculated, and C [m, j] (j = 1, 2, ..., n The optimum alignment can be obtained by tracing back from the point (m, j) that gives the minimum value of () (Non-patent Document 8).

  Here, a method is considered in which a plurality of sets of Ph ′, Mh ′, Pv ′, and Mv ′ are simultaneously calculated in bit parallel from a plurality of sets of Ph, Mh, Pv, and Mv to speed up the processing. FIG. 3 shows the positional relationship of these Boolean variables within the DP matrix. These Boolean variables represent information about the difference between the cost values C [i, j] at the lattice points 301 of the DP matrix adjacent in the horizontal direction or the vertical direction around the cell 109 in the DP matrix. From (Equation 6), the double index of the Boolean variables Ph, Mh, Pv, and Mv on the input side 302 is [i−1, j] or [i, j−1], and the sum is i + j−1. The double index of the Boolean variables Ph ′, Mh ′, Pv ′, Mv ′ on the output side 303 is [i, j] and the sum is i + j. Therefore, the values of all the Boolean variables Ph, Mh, Pv, and Mv can be calculated by calculating in order from the smallest sum k of the double indexes [i, j]. Further, it will be described below that there is a method in which these calculations can be performed in parallel without any restriction on the query length.

  By the way, the number of double indexes [i, j] whose sum is equal to k increases with k, but the double index of the Boolean variable to be referred to when tracing back in calculating the optimal alignment of similar sequences is the diagonal i = Limited to the range around j. Therefore, a method of calculating the relational expression of (Equation 7) in bit parallel is limited to this range.

Let w be the bit length of the integer word used in the basic calculation of the computer. In general, w is a power of 2, and w = 64 in the current mainstream computer. The half-word length is w ′ = w / 2. Considering the range of −w <i−j ≦ w as the range around the diagonal, the double index [i, j] whose sum is k is
If k is an even number k = 2k '
(Equation 8)
[i, j] = [k '+ w', k'−w '], ..., [k' + 1, k'−1], [k ', k'], [k'−1, k '+ 1],
..., [k'−w '+ 1, k' + w'−1]
If k is odd k = 2k '+ 1,
(Equation 9)
[i, j] = [k '+ w', k'−w '+ 1], ..., [k' + 1, k '], [k', k '+ 1],
..., [k'−w '+ 1, k' + w ']
It becomes. Therefore, if the (single) index α = k′−i is made to correspond to these double indexes (grid points), regardless of whether k is an even number or an odd number,
α = −w ', −w' + 1, ..., −1, 0, 1, ..., w'−1
Corresponding to Therefore, each bit of the w-bit word variable is indexed with these in order from the lower order bit to correspond. In the following, for simplicity of explanation, w = 4 is used, and the corresponding relationship will be specifically described with reference to the drawings.

  Figure 4 shows the correspondence between the double index [i, j] and the index of the 4-bit word variable with the sum of k in the peripheral area of the diagonal line -w <i-j ≤ w when w = 4 Is illustrated. 101 is a square lattice that represents the entire double index [i, j] of the DP matrix. The indexes i and j displayed in 102 and 103 are the indexes (character positions) of the query sequence and the target sequence, respectively. The peripheral region −w <i−j ≦ w of the diagonal line i = j shown by 104 is a range 106 sandwiched between two broken lines shown by 105. Within this range 106, the one where the sum of the double indexes is k = 5 (odd number) is at the position of four white circles (408) in the area surrounded by the broken line 407. Similarly, in this range 106, the one where the sum of the double indexes is k = 6 (even number) is at the position of four black circles (410) in the area surrounded by the broken line 409. The same applies to k = 7 (odd number). The corresponding index α = −2, −1, 0, 1 is attached to the side. The bit positions of 4-bit integers corresponding to these indexes are shown in 411.

  Similarly, the (single) index α of the word variable is made to correspond to the Boolean variables Ph, Mh, Pv, and Mv having the double index [i, j] whose sum is k. The Boolean variables Ph and Mh represent the information of the difference value of the cost of the grid points adjacent in the horizontal direction, and correspond to the index α of the left grid point (see 501 in FIG. 5). That is, Ph [i, j], Mh [i, j] is represented as Ph [α], Mh [α] by associating an index α corresponding to [i, j−1]. Similarly, the Boolean variables Pv and Mv represent the information of the difference value of the cost of the lattice points adjacent in the vertical direction, and correspond to the index α of the lower lattice point (see 502 in FIG. 5). That is, Pv [i, j] and Mv [i, j] are represented as Pv [α] and Mv [α] by associating an index α corresponding to [i, j]. FIG. 5 specifically shows an example in the case of w = 4. h [α] represents a Boolean variable Ph [α] or Mh [α], and v [α] represents a Boolean variable Pv [α] or Mv [α].

  Since there are w Ph [α], Mh [α], Pv [α], and Mv [α] for a certain k in the diagonal peripheral area 106, the bit value 503 is set to 1 according to the corresponding index α. The word variables 504 and 505 that can be packed into a word and are thus obtained are denoted as Ph, Mh, Pv, and Mv. FIG. 5 shows the case where w = 4. Here, h represents Ph or Mh, and v represents Pv or Mv. Similarly, for the Boolean variable Eq = Eq [i, j], the index α corresponding to the double index [i, j] is represented as Eq [α], and the bit value is packed into the word variable accordingly. Write things as Eq.

  These word variables can be computed in ascending order for k. FIG. 1 shows the dependency of word variable data in the calculation process when w = 4 and the position of the Boolean variable calculated at that time in the DP matrix. In the diagonal surrounding area 106 (−w <i−j ≦ w), when k is an odd number, the horizontal corresponding to the Boolean variables Ph [α], Mh [α], Pv [α], Mv [α] A side connecting adjacent lattice points in the vertical direction is represented by a thick solid line 107, and when k is an even number, a corresponding similar side is represented by a thick dotted line 108. Further, each cell 109 surrounded by the edge corresponds to the Boolean variable Eq [α]. These Boolean variables packed into word variables are Ph, Mh, Pv, Mv, and Eq, and these constitute the word data shown at 110. These are calculated in ascending order for k by the bit parallel operation 111 described below.

In order to explain the bit parallel operation 111, returning to FIG. 5, the formula to be calculated around each cell 109 is confirmed. However, the situation differs depending on whether the value of k corresponding to the left side and the upper side of the cell is even or odd. First, when k is an even number, as shown on the left of FIG. 6, the input-side Boolean variables 601 on the left and upper sides corresponding to k are Pv [α], Mv [α], Ph [α], Mh [α] And the output Boolean variables 602 on the lower and right sides corresponding to k + 1 are Ph [α], Mh [α], Pv [α], Mv [α], all of which have the same index α. Yes. Therefore, in this case, for the word variables Ph, Mh, Pv, Mv, those corresponding to k are Ph, Mh, Pv, Mv, and those corresponding to k + 1 are Ph ', Mh', Pv. If it is distinguished by writing ', Mv', the same relational expression as (Equation 7),
(Equation 4)
Xv = Eq | Mv; Pv '= Mh | ~ (Xv | Ph); Mv' = Ph &Xv;
Xh = Eq | Mh; Ph '= Mv | ~ (Xh | Pv); Mh' = Pv &Xh;
Can be calculated in bit parallel. Here, for word variables a and b, a | b represents a bitwise OR operation, a & b represents a bitwise AND operation, and ~ a represents a bitwise NOT operation. On the other hand, when k is an odd number, as shown on the right side of FIG. 6, the input side Boolean variables 603 corresponding to k are Pv [α], Mv [α], Ph [α + 1], Mh [ α + 1], and the output side Boolean variables 604 corresponding to k + 1 are Ph [α], Mh [α], Pv [α + 1], Mv [α + 1] The index α is shifted by +1. Therefore, in this case, if a bit shift operation is performed before and after the above word variables Ph, Mh, Pv, Mv and Ph ', Mh', Pv ', Mv', the same relationship as in (Equation 7) is obtained. The formula becomes applicable. That is, in this case,
(Equation 5)
Ph >> = 1; Mh >> = 1;
Xv = Eq | Mv; Pv '= Mh | ~ (Xv | Ph); Mv' = Ph &Xv;
Xh = Eq | Mh; Ph '= Mv | ~ (Xh | Pv); Mh' = Pv &Xh;
Pv '<< = 1; Mv'<< = 1;
Can be calculated in bit parallel. Here, for the word variable a, a << = 1 represents a 1-bit left shift operation, and a >> = 1 represents a 1-bit right shift operation.

On the other hand, for the word variable Eq,
If k is an even number k = 2k '
(Equation 10)
Eq (2k ') = [Eq [−w'], ..., Eq [−1], Eq [0], Eq [1], ..., Eq [w'−1]]
= [Eq [k'−w '+ 1, k' + w'−1], ...,
Eq [k'-1, k '+ 1], Eq [k', k '], Eq [k' + 1, k'-1],
..., Eq [k '+ w', k'−w ']]
If k is odd k = 2k '+ 1,
(Equation 11)
Eq (2k '+ 1) = [Eq [−w'], ..., Eq [−1], Eq [0], Eq [1], ..., Eq [w'−1]]
= [Eq [k'−w '+ 1, k' + w '], ...,
Eq [k ', k' + 1], Eq [k '+ 1, k'],
..., Eq [k '+ w', k'−w '+ 1]]
It becomes. Here, the outer [,, ...,] represents an integer value formed by arranging bits (the right end is the least significant bit), and Eq [i, j] is the i of the query array q = q1 q2 ... qm This represents a Boolean variable indicating whether or not the j th character (base) tj of the th character (base) qi and the target sequence t = t1 t2. The low-order bits of the four types of bases A, C, G, and T, Aa00, Ca01, Ga10, and Ta11 are represented as l and the high-order bits as u, and l (A) = l (G) = 0, l ( G) = l (T) = 1, u (A) = u (C) = 0, u (G) = u (T) = 1,
(Equation 12)
Eq [i, j] = not ((l (qi) xor l (tj)) or (u (qi) xor u (tj)))
Holds. Here, xor represents an exclusive OR of Boolean variables. Therefore, Tl, Tu, Boolean variables l (qi), u (qi) can be arranged in descending order from the most significant bit in Boolean variables l (tj), u (tj) in ascending order from the least significant bit. If the variables are Ql and Qu,
The word variable Eq is
(Equation 13)
Eq = ~ ((Ql ^ Tl) | (Qu ^ Tu))
It can be calculated in bit parallel like Here, for the word variables a and b, a ^ b represents a bitwise exclusive OR operation. On the other hand, Tl, Tu, Ql, and Qu can be calculated recursively for k using a bit shift operation. However, the encoded bit of the character corresponding to the query array or the target array is set to the bit vacated by the bit shift. FIG. 7 shows a calculation method of Tl, Tu, Ql, and Qu corresponding to k = 4, 5, and 6 when w = 4. In the general case, the method of updating the values of Tl, Tu, Ql, Qu when k increases is specifically expressed as an equation:
When k is odd (k = 2k '+ 1)
(Equation 14)
Ql << = 1; Ql | = l (qi);
Qu << = 1; Qu | = u (qi);
When k is an even number (k = 2k ')
(Equation 15)
Tl >> = 1; Tl | = (l (tj) <<(w−1));
Tu >> = 1; Tu | = (u (tj) <<(w−1));
It becomes. Here, qi and tj represent characters newly read from the query sequence and the target sequence, i = w ′ + k ′ + 1, j = w ′ + k ′. For word variables a and b, a = a | b is abbreviated as a | = b.

  In the above calculation, the Boolean variables Pv [α], Mv [α], Ph [α], and Mh [α] corresponding to the side in contact with the boundary 105 (i−j = ± w) of the diagonal peripheral region 106 from the inside are When calculating, it is necessary to refer to the values of the Boolean variables Pv [α], Mv [α], Ph [α], and Mh [α] corresponding to the edges that are in contact with the boundary from the outside. Since the latter values are not calculated, it is necessary to determine appropriate values in advance as boundary conditions. Therefore, it is assumed that the character comparison result between the query sequence and the target sequence is always inconsistent outside the diagonal peripheral region 106. This is a (pessimistic) assumption to avoid seeking false similar sequences. This means that C [i, j] −C [i−1, j−1] = 1 always holds outside the diagonal peripheral region 106. For this purpose, as shown in Fig. 8, it is assumed that Pv [α] = 1, Mv [α] = Ph [α] = Mh [α] = Eq [α] = 0 always holds outside the diagonal region. do it. Therefore, this is set as a boundary condition. At this time, the fact that (Equation 7) holds can be confirmed directly by substituting these values.

  When the value of k is less than or equal to the word length w, the word variables Ph, Mh, Pv, Mv are Boolean variables Pv [α], Mv [corresponding to the region where i <0, j <0 outside the DP matrix. α], Ph [α], Mh [α], Eq [α] are also referred to, but for these, Pv [α] = 1, Mv [α] = Ph [α] = By assuming that Mh [α] = Eq [α] = 0 always holds, it can be avoided that a false similar sequence is obtained.

  To start the traceback process to find the optimal alignment, it is necessary to find the minimum value of the DP cost value C [m, j] (j = 1, 2, ..., n) on the lower side of the DP matrix . For similar sequences, this minimum appears near the intersection with the diagonal. Therefore, the range for searching for the minimum value is limited to the range 902 included in the diagonal peripheral region 106 on the lower side of the DP matrix. Inside the diagonal peripheral region 106, the DP cost difference value between adjacent lattice points is calculated by the above-described method, and as shown by the column of the arrow 901 in FIG. , 0] = 0, first adding along the diagonal, and then adding in the left and right directions from the intersection with the diagonal at the lower side of the DP matrix, so that the DP matrix included in the diagonal peripheral region 106 The DP cost value on the lower side 902 can be calculated and its minimum value can be obtained. Further, the traceback processing for obtaining the optimal alignment of similar sequences refers to the difference value of the DP cost between adjacent lattice points in the diagonal peripheral region (Non-Patent Document 8). The values of Ph, Mh, Pv, and Mv are stored in the computer main memory, and these values can be referred to in the traceback process. These calculations can be performed with a small calculation cost proportional to the sequence length.

[Example 2]
In the first embodiment, the range of the double index (i, j) that satisfies −w <i−j ≦ w is adopted as the diagonal peripheral region of the DP matrix. Here, w represents the word length of the computer. In this case, in order to obtain an optimal alignment by the method of Example 1, the cumulative values of the insertion length and the deletion length must be less than w. Therefore, when it is necessary to obtain an optimal alignment that includes more insertions / deletions, it is necessary to increase the width of the region around the diagonal line of the DP matrix. In the present embodiment, a method will be described in which the method of the first embodiment is modified and applied in order to increase the width of the diagonal peripheral region of the DP matrix to twice that of the first embodiment. The same method can be used with a simple analogy when the width of the diagonal region of the DP matrix is increased by more than three times.

  When the width of the diagonal peripheral region of the DP matrix is increased to twice that of the first embodiment, the range of the double index (i, j) satisfying −2 w <i−j ≦ 2 w becomes a calculation target. Therefore, this region is divided into two partial regions sandwiching the diagonal line, −2 w <i−j ≦ 0 and 0 <i−j ≦ 2 w. Since each partial region has a width w, a bit-parallel DP matrix can be calculated by a method similar to that in the first embodiment. However, at the position where these two partial areas touch each other (touch the diagonal), the Boolean variables Pv [α], Mv [α], Ph [α], Mh [α], Eq calculated in the other partial areas. The value of [α] must be used as the boundary condition. Others can be carried out in exactly the same manner as in the first embodiment.

101 Square lattice representing the entire double index [i, j] of DP matrix
102 Query sequence index (character position)
103 Target sequence index (character position)
104 Diagonal line of DP matrix: i = j
105 Boundary line around diagonal line in DP matrix: i−j = ± w
106 Diagonal area in DP matrix: −w <i−j ≦ w
107 edge connecting adjacent grid points in horizontal or vertical direction corresponding to odd k
108 edge connecting adjacent grid points in horizontal or vertical direction corresponding to even k
109 Cell surrounded by edges connecting adjacent grid points in DP matrix
110 Word data: Ph, Mh, Pv, Mv, Eq
111 Bit parallel operation on word data
201 calculator
202 External storage device
203 Genome sequence database
204 DNA sequencer
205 DNA samples
206 Read sequence data (query sequence data)
207 Data input device
208 Parameter input device
209 Main memory
210 Seed search processor
211 Detailed comparison processing section
212 Optimal alignment calculation processor
213 Output processor
214 Output device
215 Search results (genome mapping position information)
216 Optimal alignment results (base mutation information)
301 DP matrix grid points
302 Input Boolean variable
303 Boolean variable on output side
407 Set of positions (grid points) of double index [i, j] that are included in the area surrounding the diagonal and whose sum is k = 5 (odd number)
408 Position of the double index [i, j] that is included in the area surrounding the diagonal and the sum is k = 5 (odd number) (grid point)
409 Set of positions (grid points) of double index [i, j] that are included in the area surrounding the diagonal and whose sum is k = 5 (odd number)
410 Position (grid point) of double index [i, j] that is included in the area surrounding the diagonal and whose sum is k = 5 (odd number)
411 Word variable index
501 Correspondence of (single) index α to Boolean variables corresponding to adjacent edges in the horizontal direction
502 Correspondence of (single) index α to Boolean variables corresponding to adjacent edges in the vertical direction
503 A set of Boolean variables (corresponding to adjacent edges in the horizontal or vertical direction) indexed (single) corresponding to a constant k within the region surrounding the diagonal
504 Word variable packed with a set of Boolean variables (corresponding to the adjacent edges in the horizontal direction) indexed (corresponding to the adjacent edge in the horizontal direction) corresponding to a constant k in the area surrounding the diagonal line
505 Word variable packed with a set of Boolean variables (corresponding to adjacent edges in the vertical direction) indexed (corresponding to the adjacent edge in the vertical direction) corresponding to a constant k in the area surrounding the diagonal
Boolean variable on the input side when 601 k is an even number
Boolean variable on the output side when 602 k + 1 is odd
Boolean variable on the input side when 603 k is odd
Boolean variable on the output side when 604 k + 1 is even
901 DP cost difference value accumulated to obtain DP cost value
902 Range included in the area surrounding the diagonal line 106 on the lower side of the DP matrix (the range in which the minimum DP cost value is searched)

Claims (5)

A storage device for storing a character string database;
An input device for entering a search string;
An input device for inputting an upper limit value of an alignment cost for determining that two character strings are similar;
A seed search processing unit that compares the search character string with the character string database and searches for an occurrence position in a character string database of a partial sequence shared by both (seed array);
For each occurrence position of the seed array in the character string database obtained by the seed search processing unit, a partial character string (target array) in the vicinity of the character string database and the entire search character string A detailed comparison processing unit that calculates an optimal alignment by comparison and determines whether or not a sequence similar to the entire search character string is included at the position by comparing the cost with the upper limit value;
Based on the determination result of the detailed comparison processing unit, an appearance position of the character string database where a sequence similar to the search character string appears, and an output device that outputs the optimal alignment,
The detailed comparison processing unit repeatedly performs a process (extension process) for sequentially expanding the optimal alignment for the partial sequence in order to calculate the optimal alignment for the two character string data of the search character string and the target sequence. The process based on dynamic programming,
In the decompression process, when a new character is read from each character string data and compared to decompress the optimum alignment, the absolute value of the difference in the index of the character to be compared (character position counted from the beginning of the character string) is By limiting the processing to those that are less than the word length of the computer, the number of decompression processes in which the sum of those indexes is constant is suppressed to less than the word length, and the bit parallelism of operations on the word of the computer is utilized. An optimal alignment calculation device characterized in that those decompression processes are performed in parallel by one processor. An input device for inputting a search string and a target sequence;
An input device for inputting an upper limit value of an alignment cost for determining that two character strings are similar;
By comparing the target sequence with the entire search character string to calculate an optimal alignment and comparing the cost with the upper limit value, the target sequence includes a sequence similar to the entire search character string. A detailed comparison processing unit for determining whether or not
An output device that outputs the determination result of the detailed comparison processing unit and the optimum alignment;
The detailed comparison processing unit repeatedly performs a process (extension process) for sequentially expanding the optimal alignment for the partial sequence in order to calculate the optimal alignment for the two character string data of the search character string and the target sequence. The process based on dynamic programming,
In the decompression process, when a new character is read from each character string data and compared to decompress the optimum alignment, the absolute value of the difference in the index of the character to be compared (character position counted from the beginning of the character string) is By limiting the processing to those that are less than the word length of the computer, the number of decompression processes in which the sum of those indexes is constant is suppressed to less than the word length, and the bit parallelism of operations on the word of the computer is utilized. An optimal alignment calculation device characterized in that those decompression processes are performed in parallel by one processor. In the optimal alignment calculation apparatus according to claim 1,
The cost of substitution, insertion and deletion of one character in the alignment is 1,
The optimal alignment cost between the partial character string from the beginning of the search character string to the index i and the partial character string from the beginning of the neighborhood character string to the index j is represented by C [i, j], and the equation C [i , j] −C [i, j−1] = 1, the Boolean variable is represented by Ph [i, j].
The Boolean variable Mh [i, j] represents whether the equation C [i, j] −C [i, j−1] = −1 holds,
The Boolean variable Pv [i, j] indicates whether or not the equation C [i, j] −C [i−1, j] = 1 holds.
The Boolean variable Mv [i, j] represents whether or not the equation C [i, j] −C [i−1, j] = − 1 holds.
The Boolean variable Eq (i, j) represents whether or not the character of index i of the search character string matches the character of index j of the neighboring character string,
The word length of the computer is represented by w,
For a pair of i and j whose absolute value of the difference between i and j is less than w and whose sum is k, the Boolean variable Ph [i, j] represented by 1 bit of 0 or 1 is the least significant bit for j A word that can be arranged in ascending order is represented as Ph (k), and the Boolean variables Mh [i, j], Pv [i, j], Mv [i, j], Eq [i, j] in the same manner. Are represented as Mh (k), Pv (k), Mv (k), Eq (k)
As a calculation for word variables a and b by the computer, a | b for bitwise OR operation, a & b for AND operation for each bit, ~ a for bitwise NOT operation, and a << 1 for left shift operation for 1 bit , 1-bit right shift operation is expressed as a >> = 1,
Hereinafter, for simplicity of symbols, k is fixed to one value, and the word variables Ph (k), Mh (k), Pv (k), Mv (k), Eq (k) are set to Ph, Mh. , Pv, Mv, Eq,
Xv and Xh as auxiliary word variables
When k is an even number
(Equation 1)
Xv = Eq | Mv; Pv '= Mh | ~ (Xv | Ph); Mv' = Ph &Xv;
Xh = Eq | Mh; Ph '= Mv | ~ (Xh | Pv); Mh' = Pv &Xh;
in accordance with,
When k is an odd number
(Equation 2)
Ph >> = 1; Mh >> = 1;
Xv = Eq | Mv; Pv '= Mh | ~ (Xv | Ph); Mv' = Ph &Xv;
Xh = Eq | Mh; Ph '= Mv | ~ (Xh | Pv); Mh' = Pv &Xh;
Pv '<< = 1; Mv'<< = 1;
And calculate the values of word variables Ph (k + 1), Mh (k + 1), Pv (k + 1), Mv (k + 1) by calculating Ph ', Mh', Pv ', Mv' An optimal alignment calculation device characterized in that the extension processing is calculated as described above. In the optimal alignment calculation apparatus according to claim 2,
The cost of substitution, insertion and deletion of one character in the alignment is 1,
The optimal alignment cost between the partial character string from the beginning of the search character string to the index i and the partial character string from the beginning of the neighborhood character string to the index j is represented by C [i, j], and the equation C [i , j] −C [i, j−1] = 1, the Boolean variable is represented by Ph [i, j].
The Boolean variable Mh [i, j] represents whether the equation C [i, j] −C [i, j−1] = −1 holds,
The Boolean variable Pv [i, j] indicates whether or not the equation C [i, j] −C [i−1, j] = 1 holds.
The Boolean variable Mv [i, j] represents whether or not the equation C [i, j] −C [i−1, j] = − 1 holds.
The Boolean variable Eq (i, j) represents whether or not the character of index i of the search character string matches the character of index j of the neighboring character string,
The word length of the computer is represented by w,
For a pair of i and j whose absolute value of the difference between i and j is less than w and whose sum is k, the Boolean variable Ph [i, j] represented by 1 bit of 0 or 1 is the least significant bit for j A word that can be arranged in ascending order is represented as Ph (k), and the Boolean variables Mh [i, j], Pv [i, j], Mv [i, j], Eq [i, j] in the same manner. Are represented as Mh (k), Pv (k), Mv (k), Eq (k)
As a calculation for word variables a and b by the computer, a | b for bitwise OR operation, a & b for AND operation for each bit, ~ a for bitwise NOT operation, and a << 1 for left shift operation for 1 bit , 1-bit right shift operation is expressed as a >> = 1,
Hereinafter, for simplicity of symbols, k is fixed to one value, and the word variables Ph (k), Mh (k), Pv (k), Mv (k), Eq (k) are set to Ph, Mh. , Pv, Mv, Eq,
Xv and Xh as auxiliary word variables
When k is an even number
(Equation 1)
Xv = Eq | Mv; Pv '= Mh | ~ (Xv | Ph); Mv' = Ph &Xv;
Xh = Eq | Mh; Ph '= Mv | ~ (Xh | Pv); Mh' = Pv &Xh;
in accordance with,
When k is an odd number
(Equation 2)
Ph >> = 1; Mh >> = 1;
Xv = Eq | Mv; Pv '= Mh | ~ (Xv | Ph); Mv' = Ph &Xv;
Xh = Eq | Mh; Ph '= Mv | ~ (Xh | Pv); Mh' = Pv &Xh;
Pv '<< = 1; Mv'<< = 1;
And calculate the values of word variables Ph (k + 1), Mh (k + 1), Pv (k + 1), Mv (k + 1) by calculating Ph ', Mh', Pv ', Mv' An optimal alignment calculation device characterized in that the extension processing is calculated as described above. On the computer,
Processing to accept input of search string and target sequence;
A process of accepting an input of an upper limit value of an alignment cost for determining that two character strings are similar;
By comparing the target sequence with the entire search character string to calculate an optimal alignment and comparing the cost with the upper limit value, the target sequence includes a sequence similar to the entire search character string. A detailed comparison process for determining whether or not
A program for executing the determination result of the detailed comparison process and a process of outputting the optimum alignment to an output device;
The detailed comparison process dynamically repeats a process (extension process) for sequentially expanding the optimal alignment for the partial sequence in order to calculate the optimal alignment for the two character string data of the search character string and the target sequence. Execute processing based on the programming method,
In the decompression process, when a new character is read from each character string data and compared to decompress the optimum alignment, the absolute value of the difference in the index of the character to be compared (character position counted from the beginning of the character string) is By limiting the processing to those that are less than the word length of the computer, the number of decompression processes in which the sum of those indexes is constant is suppressed to less than the word length, and the bit parallelism of operations on the word of the computer is utilized. A program characterized by performing the decompression processing in parallel by one processor.

Images