The following is a summary from the top of my head so please correct me if I'm wrong. See the Wikipedia article for more.
This is about determining the (score / edit distance of the) most similar substrings of two given strings.
For example, given the two sequences
TGTTACGG and GGTTGACTA
the optimal (w.r.t. to a certain scoring function) local alignment would be
I am using a slightly different way to score introduction and continuation of gaps than in the Wikipedia article. Technically, in addition to one scoring matrix / edit graph, you have two more - each accounting for gaps in one of the input words.
Given two input words
a = a1 a2, ..., an,
b = b1 b2 ... bm, let
s_ij denote the score of the optimal local alignment of
b1..bj. This then gives the following recurrence:
rho penalty for introducing a gap;
sigma penalty for continuing a gap; and
delta a function that assigns a score to two matched characters (e.g. 1 if equal, -1 if unequal). The higher the score the better.
This can be solved using Dynamic Programming.
On the basis of this blog article I wrote the following Haskell code:
import qualified Data.Array as Array import Data.Array scoreLocal :: String -> String -> Float scoreLocal a b = maximum mid where (m, n) = (length a, length b) -- use input as arrays for fast direct indexing v = Array.listArray (1, m) a w = Array.listArray (1, n) b -- declare scoring fns for each "edit graph level" scmid i 0 = 0 scmid 0 j = 0 scmid i j = maximum [ 0 -- start here ("free" edges from (0,0)) ,(mid ! (i-1, j-1)) + δ (v!i) (w!j) -- match or mismatch , low ! (i , j) -- end deletion / gap , upp ! (i , j) -- end insertion / gap ] scupp i 0 = 0 scupp 0 j = 0 scupp i j = maximum [ (upp ! (i, j-1)) - σ -- cont gap in v , (mid ! (i, j-1)) - ρ -- start gap in v ] sclow i 0 = 0 sclow 0 j = 0 sclow i j = maximum [ (low ! (i-1, j)) - σ -- cont gap in w , (mid ! (i-1, j)) - ρ -- start gap in w ] -- declare content of "edit graph levels" mid = Array.listArray bounds [scmid i j | (i, j) <- Array.range bounds] low = Array.listArray bounds [sclow i j | (i, j) <- Array.range bounds] upp = Array.listArray bounds [scupp i j | (i, j) <- Array.range bounds] bounds = ((0,0), (m,n)) δ v w | v == w = 1 -- match | otherwise = -2 -- mismatch ρ = 5 σ = 0.5
I dont have much experience yet with Haskell. Can (parts of this) be written in a nicer way?