trjhtr

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

GTTA-C
GTTGAC

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 a1..ai and b1..bj. This then gives the following recurrence:

With 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

Question

I dont have much experience yet with Haskell. Can (parts of this) be written in a nicer way?