align (empty) → 0.1.0.0
raw patch · 4 files changed
+204/−0 lines, 4 filesdep +basedep +uglymemodep +vectorsetup-changed
Dependencies added: base, uglymemo, vector
Files
- LICENSE +30/−0
- Setup.hs +2/−0
- align.cabal +21/−0
- src/Data/Align.hs +151/−0
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) 2014, Robin Palotai++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++ * Redistributions in binary form must reproduce the above+ copyright notice, this list of conditions and the following+ disclaimer in the documentation and/or other materials provided+ with the distribution.++ * Neither the name of Robin Palotai nor the names of other+ contributors may be used to endorse or promote products derived+ from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ align.cabal view
@@ -0,0 +1,21 @@+name: align+version: 0.1.0.0+synopsis: Sequence alignment algorithms.+description: Global or local sequence alignment, not exclusively for text.+license: BSD3+license-file: LICENSE+author: Robin Palotai+maintainer: palotai.robin@gmail.com+category: Data+build-type: Simple+cabal-version: >=1.10+source-repository head+ type: git+ location: http://bitbucket.org/robin_p/align.git++library+ exposed-modules: Data.Align+ build-depends: base >=4.6 && <5, vector >=0.10, uglymemo >= 0.1+ hs-source-dirs: src+ default-language: Haskell98+ ghc-options: -Wall
+ src/Data/Align.hs view
@@ -0,0 +1,151 @@+{-# LANGUAGE RecordWildCards #-}+module Data.Align+ ( align+ , AlignConfig+ , alignConfig+ , Step+ , Trace, traceScore, trace+ , windowedAlign+ , debugAlign+ ) where++import Data.Function (fix)+import qualified Data.List as L+import Data.MemoUgly+import Data.Ord+import qualified Data.Vector as V+import qualified Data.Vector.Generic as G+import qualified Debug.Trace as D++data AlignConfig a s = AlignConfig+ { acPairScore :: a -> a -> s+ , ac_initial_gap_penalty :: s+ , ac_gap_penalty :: s+ }++alignConfig :: (a -> a -> s) -- ^ Scoring function.+ -> s -- ^ Initial gap penalty.+ -> s -- ^ Gap penalty.+ -> AlignConfig a s+alignConfig = AlignConfig++-- | Either an unmatched item or a match.+type Step a = Either (Either a a) (a, a)++stepLeft = Left . Left+stepRight = Left . Right+stepBoth a b = Right (a,b)++isMatch :: Step a -> Bool+isMatch (Right _) = True+isMatch _ = False++isLeft :: Step a -> Bool+isLeft (Left (Left _)) = True+isLeft _ = False++isRight :: Step a -> Bool+isRight (Left (Right _)) = True+isRight _ = False++-- | The result of the alignment.+data Trace a s = Trace+ { traceScore :: s+ , trace :: [Step a]+ }++instance (Show a, Show s) => Show (Trace a s) where+ show (Trace s t) = "Trace(score = " ++ show s ++ ", steps = " ++ show t ++ ")"++(Trace s ts) `tappend` (Trace z (t:_)) = Trace (s+z) (t:ts)++-- | Utility for displaying a Char-based alignment.+debugAlign :: [Step Char] -> String+debugAlign = go [] []+ where+ go as bs [] = reverse as ++ "\n" ++ reverse bs+ go as bs (t:ts) = case t of+ Left (Left c) -> go (c:as) ('-':bs) ts+ Left (Right c) -> go ('-':as) (c:bs) ts+ Right (c, d) -> go (c:as) (d:bs) ts++-- | Aligns two sequences.+--+-- >>> :{+-- let tr = align+-- (alignConfig (\a b -> if a == b then 1 else (-0.25 :: Double)) +-- (-0.5) (-1))+-- (Data.Vector.fromList "dopple")+-- (Data.Vector.fromList "applied")+-- in do+-- print $ traceScore tr+-- putStrLn . debugAlign . trace $ tr+-- :}+-- 1.25+-- doppl-e-+-- -applied+align :: (G.Vector v a, Num s, Eq s, Ord s)+ => AlignConfig a s+ -> v a -- ^ Left sequence.+ -> v a -- ^ Right sequence.+ -> Trace a s+align (AlignConfig {..}) as bs =+ revTrace . fix (memo . go) $ (lastIndex as, lastIndex bs)+ where+ revTrace (Trace s t) = Trace s (reverse t)+ lastIndex v = G.length v - 1+ --+ go k (i,j)+ | i == (-1) || j == (-1) = + if i == j then Trace 0 []+ else if i == (-1)+ then skipInit j stepRight bs+ else skipInit i stepLeft as+ | otherwise =+ let a = as G.! i+ b = bs G.! j+ diag = k (i-1,j-1) `tappend` Trace (acPairScore a b) [stepBoth a b]+ a_gap = k (i-1, j) `tappend` Trace ac_gap_penalty [stepLeft a]+ b_gap = k ( i,j-1) `tappend` Trace ac_gap_penalty [stepRight b]+ in L.maximumBy (comparing traceScore) [diag, a_gap, b_gap]+ --+ skipInit idx stepFun xs =+ let score = ac_initial_gap_penalty * fromIntegral (idx+1)+ tr = reverse [stepFun (xs G.! xi) | xi <- [0..idx]]+ in Trace score tr++-- | Aligns long streams by performing alignment on windowed sections.+windowedAlign :: (Num s, Eq s, Ord s)+ => AlignConfig a s + -> Int -- ^ Window size.+ -> [a] -- ^ Left stream.+ -> [a] -- ^ Right stream.+ -> [Step a] -- ^ Alignment result.+windowedAlign cfg w as bs = go as bs []+ where+ go as [] rs = concat . reverse $ (map stepLeft as):rs+ go [] bs rs = concat . reverse $ (map stepRight bs):rs+ go as bs rs =+ let ahs = take w as+ bhs = take w bs+ tr0 = trace $ align cfg (V.fromList ahs) (V.fromList bhs)+ tr = let matched = dropTrailingMismatches tr0+ len = length matched+ in if len <= 1+ then matched+ -- Give opportunity for recovery given new inputs in next go.+ -- TODO this needs some experimenting with (or even option)+ else take (len `div` 2) matched+ in if null tr+ then let (a:ax) = as+ (b:bx) = bs+ -- Degenerate case, mimic some progress to avoid getting stuck.+ in go ax bx ([stepLeft a, stepRight b]:rs)+ else let ac = countMatchOr isLeft tr+ bc = countMatchOr isRight tr+ in go (drop ac as) (drop bc bs) (tr:rs)+ --+ dropTrailingMismatches =+ reverse . dropWhile (not . isMatch) . reverse+ --+ countMatchOr f = length . filter (\s -> isMatch s || f s)