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?

Improve this question
\$\endgroup\$
3
  • \$\begingroup\$ Why does the scoring function prefer the A before the gap and not after the gap? \$\endgroup\$
    – Vogel612
    May 23, 2018 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\$
    – Dan Oberlam
    May 23, 2018 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, 2018 at 9:22

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.