align-affine (empty) → 0.1.0.0
raw patch · 5 files changed
+345/−0 lines, 5 filesdep +basedep +containersdep +transformerssetup-changed
Dependencies added: base, containers, transformers, vector
Files
- CHANGELOG.md +3/−0
- LICENSE.md +28/−0
- Setup.hs +2/−0
- align-affine.cabal +25/−0
- src/Data/Align/Affine.hs +287/−0
+ CHANGELOG.md view
@@ -0,0 +1,3 @@+# [0.1.0.0] - 2021-11-26++- First version.
+ LICENSE.md view
@@ -0,0 +1,28 @@+Copyright 2021, Charles Li. Forked from copyright 2014, Robin Palotai, Google Inc.++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 the copyright holder 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-affine.cabal view
@@ -0,0 +1,25 @@+name: align-affine+version: 0.1.0.0+synopsis: Sequence alignment with an affine gap penalty model+description: Global & multi-sequence alignment using the Needleman-Wunsch & Center Star algorithms, with an affine gap penalty model.+license: BSD3+license-file: LICENSE.md+author: Charles Li+maintainer: charles.li@ehess.fr+category: Data+extra-source-files: CHANGELOG.md+build-type: Simple+cabal-version: >=1.10+source-repository head+ type: git+ location: https://github.com/chchch/align-affine.git++library+ exposed-modules: Data.Align.Affine+ build-depends: base >=4.6 && <5,+ containers >=0.6 && <0.7,+ transformers >=0.5 && <0.7,+ vector >=0.10 && <0.13+ hs-source-dirs: src+ default-language: Haskell98+ ghc-options: -Wall
+ src/Data/Align/Affine.hs view
@@ -0,0 +1,287 @@+{-# LANGUAGE RecordWildCards #-}+-- | Pairwise & multi-sequence alignment with an affine gap penalty model. Forked from "Data.Align". This algorithm uses a gap opening penalty and a gap extension penalty to score alignments.+module Data.Align.Affine+ (+ -- * Global alignment+ align+ , AlignConfig+ , alignConfig+ , Step+ , Trace, traceScore, trace+ , debugAlign, debugStrAlign+ -- * Multi-sequence alignment+ , centerStar+ , MultiStep, center, others, stepOfAll+ , MultiTrace, centerIndex, otherIndices, allIndices, multiTrace+ , debugMultiAlign+ ) where++import Control.Monad.Trans.State.Strict (evalState, gets, modify)+import Data.Function (on)+import Data.Maybe (fromMaybe)+import qualified Data.List as L+import qualified Data.Map.Strict as M+import Data.Ord (comparing)+import qualified Data.Vector.Generic as G++-- | Either an unmatched item or a match.+type Step a = Either (Either a a) (a, a)++stepLeft :: a -> Either (Either a b1) b2+stepLeft = Left . Left++stepRight :: b1 -> Either (Either a b1) b2+stepRight = Left . Right++stepBoth :: a1 -> b -> Either a2 (a1, b)+stepBoth a b = Right (a,b)++-- | 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 ++ ")"++tappend :: Num s => Trace a s -> (s, Step a) -> Trace a s+Trace a b `tappend` (y,z) = Trace (a+y) (z:b)++data AffineTrace a s = AffineTrace {+ at_max :: Trace a s+ , at_left_gap :: Trace a s+ , at_right_gap :: Trace a s+ }++data AlignConfig a s = AlignConfig+ { acPairScore :: a -> a -> s+ , acGapOpen :: s+ , acGapExtension :: s+ }++-- | This algorithm uses an affine gap penalty model. See section 12.6 "Convex gap weights" in Gusfield 1997: <https://doi.org/10.1017/CBO9780511574931>.+alignConfig :: (a -> a -> s) -- ^ Scoring function.+ -> s -- ^ Gap opening score.+ -> s -- ^ Gap extension score.+ -> AlignConfig a s+alignConfig = AlignConfig++-- | Aligns two sequences using the Needleman-Wunsch algorithm. See Needleman & Wunsch 1970: <https://doi.org/10.1016/0022-2836(70)90057-4>.+align :: (G.Vector v a, Num s, Ord s)+ => AlignConfig a s+ -> v a -- ^ Left sequence.+ -> v a -- ^ Right sequence.+ -> Trace a s+align AlignConfig{..} as bs =+ let p = (lastIndex as, lastIndex bs)+ in revTrace . at_max $ evalState (go p) M.empty+ where+ revTrace (Trace s t) = Trace s (reverse t)+ lastIndex v = G.length v - 1+ --+ go p = do+ res <- gets $ M.lookup p+ case res of+ Just r -> return r+ Nothing -> do+ newRes <- pgo p+ modify (M.insert p newRes)+ return newRes+ --+ pgo (i,j)+ | i == (-1) || j == (-1) = return $+ if i == j then AffineTrace (Trace 0 []) (Trace 0 []) (Trace 0 [])+ else if i == (-1)+ then skipInit j stepRight bs+ else skipInit i stepLeft as+ | otherwise = do+ let a = as G.! i+ b = bs G.! j+ diag <- go (i-1,j-1)+ let diag_max = (at_max diag) `tappend` (acPairScore a b, stepBoth a b)+ + a_gaps <- go (i-1, j)+ let a_gap1 = (at_max a_gaps) `tappend` (acGapOpen + acGapExtension, stepLeft a)+ let a_gap2 = (at_left_gap) a_gaps `tappend` (acGapExtension, stepLeft a)+ let a_gap_max = L.maximumBy (comparing traceScore) [a_gap1, a_gap2]+ + b_gaps <- go ( i,j-1)+ let b_gap1 = (at_max b_gaps) `tappend` (acGapOpen + acGapExtension, stepRight b)+ let b_gap2 = (at_right_gap b_gaps) `tappend` (acGapExtension, stepRight b)+ let b_gap_max = L.maximumBy (comparing traceScore) [b_gap1, b_gap2]+ + let maxi = L.maximumBy (comparing traceScore) [diag_max, a_gap_max, b_gap_max]+ return $ AffineTrace maxi a_gap_max b_gap_max+ --+ skipInit idx stepFun xs =+ let score = acGapOpen + acGapExtension * fromIntegral (idx+1)+ tr = reverse [stepFun (xs G.! xi) | xi <- [0..idx]]+ in AffineTrace (Trace score tr) (Trace score tr) (Trace score tr)++-- | Utility for displaying a Char-based alignment.+--+-- == Example+-- +-- > >>> :{+-- > let tr = align (alignConfig (\a b -> if a == b then 0.5 else (-0.5)) +-- > (-3) (-0.25)) +-- > (Data.Vector.fromList "circumambulate") +-- > (Data.Vector.fromList "perambulatory")+-- > in do+-- > print $ traceScore tr+-- > putStrLn . debugAlign . trace $ tr+-- > :}+-- +-- == Output+--+-- > -4.75+-- > circumambulate--+-- > per---ambulatory+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++-- | Utility for displaying a String-based alignment.+-- +-- == Example+--+-- >>> :{+-- let tr = align (alignConfig (\a b -> if a == b then 0.5 else (-0.5)) +-- (-3) (-0.25)) +-- (Data.Vector.fromList ["kra","ya","ṇā","ddha","ra","ṇā","tyā","cñā","yāḥ"]) +-- (Data.Vector.fromList ["bha","ra","ṇā","da","pa","ha","ra","ṇā","tyā","cña","yā"])+-- in do+-- print $ traceScore tr+-- putStrLn . debugStrAlign . trace $ tr+-- :}+--+-- == Output+--+-- > -4.0+-- > |kra|ya|ṇā|ddha|--|--|ra|ṇā|tyā|cñā|yāḥ|+-- > |bha|ra|ṇā|da |pa|ha|ra|ṇā|tyā|cñā|yā |+debugStrAlign :: [Step String] -> String+debugStrAlign = go id id+ where+ go as bs [] = as "|\n" ++ bs "|"+ go as bs (t:ts) = case t of+ Left (Left c) -> go (as . ("|"++) . (c ++)) + (bs . ("|"++) . (take (length c) (repeat '-') ++)) ts+ Left (Right c) -> go (as . ("|"++) . (take (length c) (repeat '-') ++)) + (bs . ("|"++) . (c ++)) ts+ Right (c,d) -> go (as . ("|"++) . (c ++) . (filldc ++)) + (bs . ("|"++) . (d ++) . (fillcd ++)) ts+ where+ fill n = take n $ repeat ' '+ filldc = fill (length d - length c)+ fillcd = fill (length c - length d)++-- | A step in a multi-sequence alignment.+data MultiStep a = MultiStep+ { center :: Maybe a -- ^ 'Nothing' means gap insertion.+ , others :: [Maybe a] -- ^ Parallel to 'otherIndices'.+ }++-- | The result of a multi-sequence alignment.+data MultiTrace i a s = MultiTrace+ { centerIndex :: i+ , otherIndices :: [i]+ , multiTrace :: [MultiStep a]+ }++-- | The center step followed by other steps.+stepOfAll :: MultiStep a -> [Maybe a]+stepOfAll MultiStep{..} = center:others++-- | The center index followed by other indices.+allIndices :: MultiTrace i a s -> [i]+allIndices MultiTrace{..} = centerIndex:otherIndices++-- | Align multiple sequences using the Center Star method. +--+-- See, for example, section 14.6.2, "A bounded-error approximation method for SP alignment" in Gusfield 1997.+-- +-- Assumes the list of sequences to have length > 1, and the indices to be unique.+centerStar :: (G.Vector v a, Num s, Ord s, Ord i)+ => AlignConfig a s+ -> [(i, v a)] -- TODO use internal indices rather to make uniqueness sure+ -> MultiTrace i a s+centerStar conf vs =+ let (firstPair:rest) = centerPairs+ initialTrace = MultiTrace+ { centerIndex = fst . fst $ firstPair+ , otherIndices = [snd . fst $ firstPair]+ , multiTrace = initialSteps . trace . snd $ firstPair+ }+ in foldl mergePair initialTrace rest+ where+ initialSteps = go []+ where+ go acc [] = reverse acc+ go acc (s:xs) = go (conv s []:acc) xs+ --+ conv s rest = case s of+ Right (c, d) -> MultiStep (Just c) (Just d:rest)+ Left (Left c) -> MultiStep (Just c) (Nothing:rest)+ Left (Right d) -> MultiStep Nothing (Just d:rest)+ --+ mergePair MultiTrace{..} ((_,j), tr) = MultiTrace+ { centerIndex = centerIndex+ , otherIndices = j:otherIndices+ , multiTrace = mergeSteps multiTrace (trace tr)+ }+ where+ mergeSteps mss' = go [] mss'+ where+ noOthers = map (const Nothing) . others . head $ mss'+ --+ go acc [] [] = reverse acc+ go acc (MultiStep{..}:mss) [] =+ go (MultiStep center (Nothing:others):acc) mss []+ go acc [] (s:ss) = go (conv s noOthers:acc) [] ss+ go acc (m@MultiStep{..}:mss) (s:ss) = case (center, s) of+ (Nothing, Left (Right d)) ->+ go (MultiStep center (Just d:others):acc) mss ss+ (Nothing, _) ->+ go (MultiStep center (Nothing:others):acc) mss (s:ss)+ (Just _, Right (_, d)) ->+ go (MultiStep center (Just d:others):acc) mss ss+ (Just _, Left (Left _)) ->+ go (MultiStep center (Nothing:others):acc) mss ss+ (Just _, Left (Right d)) ->+ go (MultiStep Nothing (Just d:noOthers):acc) (m:mss) ss+ --+ centerPairs+ = snd -- drop cache+ . L.maximumBy (comparing fst)+ . map (\g -> (starSum g, g)) -- cache scores+ . L.groupBy ((==) `on` (fst . fst))+ . L.sortBy (comparing fst)+ $ pairAligns+ where+ pairAligns = do+ ((i,v):rest) <- L.tails vs+ (j,w) <- rest+ let tr = align conf v w+ [((i,j), tr), ((j,i), flipLR tr)]+ where+ flipLR tr = tr { trace = map go . trace $ tr }+ where+ go (Left (Left a)) = Left (Right a)+ go (Left (Right a)) = Left (Left a)+ go (Right (c,d)) = Right (d,c) + --+ starSum = sum . map (traceScore . snd)++-- | Renders a char-based multi-alignment result to a string.+debugMultiAlign :: [MultiStep Char] -> String+debugMultiAlign =+ unlines . map (map charOrDash) . L.transpose . map stepOfAll+ where+ charOrDash = fromMaybe '-'