2
\$\begingroup\$

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:

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?

\$\endgroup\$
  • \$\begingroup\$ Why does the scoring function prefer the A before the gap and not after the gap? \$\endgroup\$ – Vogel612 May 23 '18 at 13:16
  • 2
    \$\begingroup\$ Your second question, of how to modify the code to compute that value, is off-topic; we only review existing code, we don't enhance or provide enhancements to the functionality. \$\endgroup\$ – Dannnno May 23 '18 at 13:48
  • \$\begingroup\$ I removed the off-topic part of your question, it was either that or getting the question closed. Please take a look at the help center. Feature requests are not something we do. \$\endgroup\$ – Mast May 24 '18 at 9:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.