Hirschberg's algorithm

Algorithm information

Hirschberg's algorithm is a generally applicable algorithm for optimal sequence alignment. BLAST and FASTA are suboptimal heuristics. If x and y are strings, where length(x) = n and length(y) = m, the Needleman-Wunsch algorithm finds an optimal alignment in O(nm) time, using O(nm) space. Hirschberg's algorithm is a clever modification of the Needleman-Wunsch Algorithm which still takes O(nm) time, but needs only O(min{n,m}) space. One application of the algorithm is finding sequence alignments of DNA or protein sequences. It is also a space-efficient way to calculate the longest common subsequence between two sets of data such as with the common diff tool.

The Hirschberg algorithm can be derived from the Needleman-Wunsch algorithm by observing that:

1. one can compute the optimal alignment score by only storing the current and previous row of the Needleman-Wunsch score matrix;
2. if ( Z , W ) = NW ⁡ ( X , Y ) is the optimal alignment of ( X , Y ) , and X = X l + X r is an arbitrary partition of X , there exists a partition Y l + Y r of Y such that NW ⁡ ( X , Y ) = NW ⁡ ( X l , Y l ) + NW ⁡ ( X r , Y r ) .

Algorithm description

X i denotes the i-th character of X , where 1 < i ⩽ length ⁡ ( X ) . X i : j denotes a substring of size j − i + 1 , ranging from i-th to the j-th character of X . rev ⁡ ( X ) is the reversed version of X .

X and Y are sequences to be aligned. Let x be a character from X , and y be a character from Y . We assume that Del ⁡ ( x ) , Ins ⁡ ( y ) and Sub ⁡ ( x , y ) are well defined integer-valued functions. These functions represent the cost of deleting x , inserting y , and replacing x with y , respectively.

We define NWScore ⁡ ( X , Y ) , which returns the last line of the Needleman-Wunsch score matrix S c o r e ( i , j ) :

function NWScore(X,Y) Score(0,0) = 0 for j=1 to length(Y) Score(0,j) = Score(0,j-1) + Ins(Y_j) for i=1 to length(X) Score(i,0) = Score(i-1,0) + Del(X_i) for j=1 to length(Y) scoreSub = Score(i-1,j-1) + Sub(X_i, Y_j) scoreDel = Score(i-1,j) + Del(X_i) scoreIns = Score(i,j-1) + Ins(Y_j) Score(i,j) = max(scoreSub, scoreDel, scoreIns) end end for j=0 to length(Y) LastLine(j) = Score(length(X),j) return LastLine

Note that at any point, NWScore only requires the two most recent rows of the score matrix. Thus, NWScore can be implemented in O ( min ⁡ { length ⁡ ( X ) , length ⁡ ( Y ) } ) space, although the algorithm above takes O(i*j) space.

The Hirschberg algorithm follows:

function Hirschberg(X,Y) Z = "" W = "" if length(X) == 0 for i=1 to length(Y) Z = Z + '-' W = W + Y_i end else if length(Y) == 0 for i=1 to length(X) Z = Z + X_i W = W + '-' end else if length(X) == 1 or length(Y) == 1 (Z,W) = NeedlemanWunsch(X,Y) else xlen = length(X) xmid = length(X)/2 ylen = length(Y) ScoreL = NWScore(X_1:xmid, Y) ScoreR = NWScore(rev(X_xmid+1:xlen), rev(Y)) ymid = arg max ScoreL + rev(ScoreR) (Z,W) = Hirschberg(X_1:xmid, y_1:ymid) + Hirschberg(X_xmid+1:xlen, Y_ymid+1:ylen) end return (Z,W)

In the context of Observation (2), assume that X l + X r is a partition of X . Index y m i d is computed such that Y l = Y 1 : y m i d and Y r = Y y m i d + 1 : length ⁡ ( Y ) .

Example

Let

X = A G T A C G C A , Y = T A T G C , Del ⁡ ( x ) = − 2 , Ins ⁡ ( y ) = − 2 , Sub ⁡ ( x , y ) = { + 2 , if  x = y − 1 , if  x ≠ y .

The optimal alignment is given by

W = AGTACGCA Z = --TATGC-

Indeed, this can be verified by backtracking its corresponding Needleman-Wunsch matrix:

T A T G C 0 -2 -4 -6 -8 -10 A -2 -1 0 -2 -4 -6 G -4 -3 -2 -1 0 -2 T -6 -2 -4 0 -2 -1 A -8 -4 0 -2 -1 -3 C -10 -6 -2 -1 -3 1 G -12 -8 -4 -3 1 -1 C -14 -10 -6 -5 -1 3 A -16 -12 -8 -7 -3 1

One starts with the top level call to Hirschberg ⁡ ( A G T A C G C A , T A T G C ) , which splits the first argument in half: X = A G T A + C G C A . The call to NWScore ⁡ ( A G T A , Y ) produces the following matrix:

T A T G C 0 -2 -4 -6 -8 -10 A -2 -1 0 -2 -4 -6 G -4 -3 -2 -1 0 -2 T -6 -2 -4 0 -2 -1 A -8 -4 0 -2 -1 -3

Likewise, NWScore ⁡ ( rev ⁡ ( C G C A ) , rev ⁡ ( Y ) ) generates the following matrix:

C G T A T 0 -2 -4 -6 -8 -10 A -2 -1 -3 -5 -4 -6 C -4 0 -2 -4 -6 -5 G -6 -2 2 0 -2 -4 C -8 -4 0 1 -1 -3

Their last lines (after reversing the latter) and sum of those are respectively

ScoreL = [ -8 -4 0 -2 -1 -3 ] rev(ScoreR) = [ -3 -1 1 0 -4 -8 ] Sum = [-11 -5 1 -2 -5 -11]

The maximum (shown in bold) appears at ymid = 2, producing the partition Y = T A + T G C .

The entire Hirschberg recursion (which we omit for brevity) produces the following tree:

(AGTACGCA,TATGC) / (AGTA,TA) (CGCA,TGC) / / (AG, ) (TA,TA) (CG,TG) (CA,C) / / (T,T) (A,A) (C,T) (G,G)

The leaves of the tree contain the optimal alignment.