diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -2,6 +2,15 @@ ## Unreleased changes +## [0.1.1.1] - 2020-03-24 +### Fixed +- Affine alignment matrix construction and traceback. +### Added +- `similarity'` family of functions to calculate similarity directly on `AlignnmentResult`. +### Changed +- Generalize affine gap matrix construction; +- Improve alignment documentation. + ## [0.1.1.0] - 2019-06-17 ### Added - Typeclass `IsGap`. diff --git a/cobot.cabal b/cobot.cabal --- a/cobot.cabal +++ b/cobot.cabal @@ -4,10 +4,10 @@ -- -- see: https://github.com/sol/hpack -- --- hash: 5253d8bdc8faf78e15758bbf5ca466c06c743718a111840ff74966c4932315ad +-- hash: 19558d31d73118631f947a34b8de243404cf58e8942adbd414bf445b21381b86 name: cobot -version: 0.1.1.0 +version: 0.1.1.1 synopsis: Computational biology toolkit to collaborate with researchers in constructive protein engineering description: Please see the README on GitHub at category: Bio diff --git a/src/Bio/Chain/Alignment.hs b/src/Bio/Chain/Alignment.hs --- a/src/Bio/Chain/Alignment.hs +++ b/src/Bio/Chain/Alignment.hs @@ -1,14 +1,26 @@ module Bio.Chain.Alignment - ( AlignmentResult (..), SimpleGap, SimpleGap2, AffineGap (..), AffineGap2, Operation (..) + ( + -- * Alignment function + align + -- ** Alignment algorithms , EditDistance (..) , GlobalAlignment (..), LocalAlignment (..), SemiglobalAlignment (..) + , SimpleGap, SimpleGap2, AffineGap (..), AffineGap2 , IsGap (..) - , align + -- ** Alignment result types + , AlignmentResult (..), Operation (..) + -- * Viewing alignment results , viewAlignment , prettyAlignmment + -- * Similarity functions + -- $similarity + , similarityGen' , similarityGen + , differenceGen' , differenceGen + , similarity' , similarity + , difference' , difference ) where @@ -82,69 +94,69 @@ -> Index m -> Index m' -> [Operation (Index m) (Index m')] -traceback algo mat s t i' j' = helper i' j' [] +traceback algo mat s t i' j' = helper i' j' Match [] where - helper i j ar | isStop (cond algo) mat s t i j = ar - | isVert (cond algo) mat s t i j = helper (pred i) j (DELETE (pred i):ar) - | isHoriz (cond algo) mat s t i j = helper i (pred j) (INSERT (pred j):ar) - | isDiag (cond algo) mat s t i j = helper (pred i) (pred j) (MATCH (pred i) (pred j):ar) - | otherwise = error "Alignment traceback: you cannot be here" + helper i j prevOp ar + | isStop (cond algo) mat s t i j = ar + | otherwise = + let (nextOp, nextI, nextJ, op) = doMove (cond algo) mat s t i j prevOp + in helper nextI nextJ nextOp $ op:ar ---------------------------------------------------------------------------------------------------------- - -- - -- Some TIPS for using the functions below - -- - -- These are generic variants of similarity and difference functions alongside with their specialised variants. - -- Generic versions take the alignment algorithm used for sequence alignment, - -- an equality function on elements of both sequences to calculate hamming distance on aligned sequences, - -- and the sequences themselves. - -- - -- Sample usage of generic functions: - -- - -- > similarityGen (GlobalAlignment (\x y -> if x == ord y then 1 else 0) (AffineGap (-11) (-1))) (\x y -> x == ord y) [ord 'R'.. ord 'z'] ['a'..'z'] - -- > 0.63414633 - -- - -- This one will calculate similarity between a list if `Int`s and a list of `Char`s. - -- Generic scoring function used in alignment is `\x y -> if x == ord y then 1 else 0` - -- Generic equality function used in hamming distance is `\x y -> x == ord y` - -- - -- - -- Specialised versions do not take the equality function as the sequences are already constrained to have `Eq` elements. - -- - -- Sample usage of specialised function is the same as before: - -- - -- > seq1 :: String - -- > seq1 = "EVQLLESGGGLVQPGGSLRLSCAASGFTFSSYAMSWVRQAPGKGLEWVSAISGSGGSTYYADSVKGRFTISRDNSKNTLYLQMNSLRAEDTAVYYCAKVQLERYFDYWGQGTLVTVSS" - -- > - -- > seq2 :: String - -- > seq2 = "EVQLLESGGGLVQPGGSLRLSAAASGFTFSTFSMNWVRQAPGKGLEWVSYISRTSKTIYYADSVKGRFTISRDNSKNTLYLQMNSLRAEDTAVYYVARGRFFDYWGQGTLVTVS" - -- > - -- > similarity (GlobalAlignment blosum62 (AffineGap (-11) (-1))) s1 s2 - -- > 0.8130081 - -- - -- Sometimes for biological reasons gaps appearing in one of two sequences, that are being aligned, - -- are not physical. For that reason we might want to use different gap penalties when aligning these sequences. - -- - -- Example of usage of different gaps when aligning two sequences is presented below: - -- - -- > seq1 :: String - -- > seq1 = "AAAAAAGGGGGGGGGGGGTTTTTTTTT" - -- > - -- > seq2 :: String - -- > seq2 = "AAAAAATTTTTTTTT" - -- > - -- > gapForSeq1 :: AffineGap - -- > gapForSeq1 = AffineGap (-5) (-1) - -- > - -- > gapForSeq2 :: AffineGap - -- > gapForSeq2 = AffineGap (-1000) (-1000) -- basically, we forbid gaps on @seq2@ - -- > - -- > local = LocalAlignment nuc44 (gapForSeq1, gapForSeq2) - -- > - -- > viewAlignment (align local seq1 seq2) == ("TTTTTTTTT", "TTTTTTTTT") - -- ---------------------------------------------------------------------------------------------------------- +{- $similarity +These are generic variants of similarity and difference functions alongside with their specialised variants. +Generic versions take the alignment algorithm used for sequence alignment, +an equality function on elements of both sequences to calculate hamming distance on aligned sequences, +and the sequences themselves. +Sample usage of generic functions: + +>>> similarityGen (GlobalAlignment (\x y -> if x == ord y then 1 else 0) (AffineGap (-11) (-1))) (\x y -> x == ord y) [ord 'R'.. ord 'z'] ['a'..'z'] +0.63414633 + +This one will calculate similarity between a list of 'Int's and a list of 'Char's. +Generic scoring function used in alignment is @\\x y -> if x == ord y then 1 else 0@. +Generic equality function used in hamming distance is @\\x y -> x == ord y@. + +Specialised versions do not take the equality function as the sequences are already constrained to have 'Eq' elements. + +Sample usage of specialised function is the same as before: + +>>> :{ +seq1 :: String +seq1 = "EVQLLESGGGLVQPGGSLRLSCAASGFTFSSYAMSWVRQAPGKGLEWVSAISGSGGSTYYADSVKGRFTISRDNSKNTLYLQMNSLRAEDTAVYYCAKVQLERYFDYWGQGTLVTVSS" + +seq2 :: String +seq2 = "EVQLLESGGGLVQPGGSLRLSAAASGFTFSTFSMNWVRQAPGKGLEWVSYISRTSKTIYYADSVKGRFTISRDNSKNTLYLQMNSLRAEDTAVYYVARGRFFDYWGQGTLVTVS" + +similarity (GlobalAlignment blosum62 (AffineGap (-11) (-1))) s1 s2 +:} +0.8130081 + +Sometimes for biological reasons gaps appearing in one of two sequences, that are being aligned, +are not physical. For that reason we might want to use different gap penalties when aligning these sequences. + +Example of usage of different gaps when aligning two sequences is presented below: + +>>> :{ +seq1 :: String +seq1 = "AAAAAAGGGGGGGGGGGGTTTTTTTTT" + +seq2 :: String +seq2 = "AAAAAATTTTTTTTT" + +gapForSeq1 :: AffineGap +gapForSeq1 = AffineGap (-5) (-1) + +gapForSeq2 :: AffineGap +gapForSeq2 = AffineGap (-1000) (-1000) -- basically, we forbid gaps on @seq2@ + +local = LocalAlignment nuc44 (gapForSeq1, gapForSeq2) + +viewAlignment (align local seq1 seq2) +:} +("TTTTTTTTT", "TTTTTTTTT") +-} + -- | Calculate similarity and difference between two sequences, aligning them first using given algorithm. -- similarityGen :: forall algo m m'.(SequenceAlignment algo, Alignable m, Alignable m') @@ -153,9 +165,18 @@ -> m -> m' -> R -similarityGen algo genericEq s t = fromIntegral hamming / fromIntegral len +similarityGen algo genericEq s t = similarityGen' (align algo s t) genericEq + +-- | Calculate similarity by precomputed 'AlignmentResult'. +similarityGen' :: forall m m'. (Alignable m, Alignable m') + => AlignmentResult m m' + -> (IxValue m -> IxValue m' -> Bool) + -> R +similarityGen' res genericEq = fromIntegral hamming / fromIntegral len where - operations = alignment (align algo s t) + operations = alignment res + s = sequence1 res + t = sequence2 res len = length operations hamming = sum $ toScores <$> operations @@ -170,6 +191,10 @@ -> R similarity algo = similarityGen algo (==) +similarity' :: forall m m'.(Alignable m, Alignable m', IxValue m ~ IxValue m', Eq (IxValue m), Eq (IxValue m')) + => AlignmentResult m m' + -> R +similarity' res = similarityGen' res (==) differenceGen :: forall algo m m'.(SequenceAlignment algo, Alignable m, Alignable m') => algo (IxValue m) (IxValue m') @@ -179,6 +204,11 @@ -> R differenceGen algo genericEq s t = 1.0 - similarityGen algo genericEq s t +differenceGen' :: forall m m'.(Alignable m, Alignable m') + => AlignmentResult m m' + -> (IxValue m -> IxValue m' -> Bool) + -> R +differenceGen' res genericEq = 1.0 - similarityGen' res genericEq difference :: forall algo m m'.(SequenceAlignment algo, Alignable m, Alignable m', IxValue m ~ IxValue m', Eq (IxValue m), Eq (IxValue m')) => algo (IxValue m) (IxValue m') @@ -187,7 +217,12 @@ -> R difference algo = differenceGen algo (==) --- | View alignment results as simple strings with gaps +difference' :: forall m m'.(Alignable m, Alignable m', IxValue m ~ IxValue m', Eq (IxValue m), Eq (IxValue m')) + => AlignmentResult m m' + -> R +difference' res = differenceGen' res (==) + +-- | View alignment results as simple strings with gaps. -- viewAlignment :: forall m m'.(Alignable m, Alignable m', Symbol (IxValue m), Symbol (IxValue m')) => AlignmentResult m m' -> (String, String) viewAlignment ar = unzip (toChars <$> alignment ar) diff --git a/src/Bio/Chain/Alignment/Algorithms.hs b/src/Bio/Chain/Alignment/Algorithms.hs --- a/src/Bio/Chain/Alignment/Algorithms.hs +++ b/src/Bio/Chain/Alignment/Algorithms.hs @@ -18,8 +18,8 @@ import Data.Array.Unboxed (Ix (..), UArray, (!)) --- | Alignnment methods --- +-- Alignnment methods + newtype EditDistance e1 e2 = EditDistance (e1 -> e2 -> Bool) data GlobalAlignment a e1 e2 = GlobalAlignment (Scoring e1 e2) a data LocalAlignment a e1 e2 = LocalAlignment (Scoring e1 e2) a @@ -57,23 +57,57 @@ (lowerT, _) = bounds t in i == lowerS || j == lowerT || m' ! (i, j, Match) == 0 -horiz :: (Alignable m, Alignable m', IsGap g) => g -> Matrix m m' -> m -> m' -> Index m -> Index m' -> Bool -horiz g m s t i j = (j > lowerT) && ((i == lowerS) || (m ! (i, pred j, Match) + add == m ! (i, j, Match))) - where - add | isAffine g = m ! (i, pred j, Insert) - | otherwise = insertCostOpen g - - (lowerT, _) = bounds t - (lowerS, _) = bounds s +{-# INLINE move #-} +move + :: (Alignable m, Alignable m', IsGap g) + => g + -> Move m m' +move g + | isAffine g = moveAffine g + | otherwise = moveSimple g -vert :: (Alignable m, Alignable m', IsGap g) => g -> Matrix m m' -> m -> m' -> Index m -> Index m' -> Bool -vert g m s t i j = (i > lowerS) && ((lowerT == j) || (m ! (pred i, j, Match) + add == m ! (i, j, Match))) - where - add | isAffine g = m ! (pred i, j, Delete) - | otherwise = deleteCostOpen g +{-# INLINE moveSimple #-} +moveSimple + :: (Alignable m, Alignable m', IsGap g) + => g + -> Move m m' +moveSimple g m s t i j _ + | i == lowerS = (Match, lowerS, pred j, INSERT $ pred j) + | j == lowerT = (Match, pred i, lowerT, DELETE $ pred i) + | m ! (pred i, j, Match) + deleteCostOpen g == m ! (i, j, Match) = (Match, pred i, j, DELETE $ pred i) + | m ! (i, pred j, Match) + insertCostOpen g == m ! (i, j, Match) = (Match, i, pred j, INSERT $ pred j) + | otherwise = (Match, pred i, pred j, MATCH (pred i) (pred j)) + where + (lowerT, _) = bounds t + (lowerS, _) = bounds s - (lowerT, _) = bounds t - (lowerS, _) = bounds s +{-# INLINE moveAffine #-} +-- | Move function for affine alignment traceback. Implements a "Manhattan grid". +-- +-- See here: +-- or file @doc/Affine_Alignment.pdf@ in the repository. +moveAffine + :: (Alignable m, Alignable m', IsGap g) + => g + -> Move m m' +moveAffine g m s t i j prevOp + | i == lowerS = (Insert, lowerS, pred j, INSERT $ pred j) + | j == lowerT = (Delete, pred i, lowerT, DELETE $ pred i) + | otherwise = + case prevOp of + Delete + | m ! (i, j, Delete) == m ! (pred i, j, Delete) + deleteCostExtend g -> (Delete, pred i, j, DELETE $ pred i) + | otherwise -> (Match, pred i, j, DELETE $ pred i) + Insert + | m ! (i, j, Insert) == m ! (i, pred j, Insert) + insertCostExtend g -> (Insert, i, pred j, INSERT $ pred j) + | otherwise -> (Match, i, pred j, INSERT $ pred j) + Match + | m ! (i, j, Match) == m ! (i, j, Delete) -> move g m s t i j Delete + | m ! (i, j, Match) == m ! (i, j, Insert) -> move g m s t i j Insert + | otherwise -> (Match, pred i, pred j, MATCH (pred i) (pred j)) + where + (lowerT, _) = bounds t + (lowerS, _) = bounds s -- | Default condition of moving diagonally in traceback. -- @@ -141,7 +175,7 @@ -- Conditions of traceback are described below -- {-# INLINE cond #-} - cond (GlobalAlignment subC gap) = Conditions defStop (defDiag subC) (vert gap) (horiz gap) + cond (GlobalAlignment _ gap) = Conditions defStop (move gap) -- Start from bottom right corner -- @@ -199,22 +233,7 @@ writeArray matrix (ixS, ixT, Match) $ deleteCostOpen g + (deleteCostExtend g) * pred (index (lowerS, nilS) ixS) writeArray matrix (ixS, ixT, Delete) $ deleteCostExtend g writeArray matrix (ixS, ixT, Insert) $ insertCostOpen g - | otherwise -> do - predDiag <- matrix `readArray` (pred ixS, pred ixT, Match) - predS <- matrix `readArray` (pred ixS, ixT, Match) - predT <- matrix `readArray` ( ixS, pred ixT, Match) - - delCost <- matrix `readArray` (pred ixS, ixT, Delete) - insCost <- matrix `readArray` ( ixS, pred ixT, Insert) - - let maxScore = maximum [ predDiag + sub ixS ixT - , predS + delCost - , predT + insCost - ] - - writeArray matrix (ixS, ixT, Delete) $ if predS + delCost == maxScore then deleteCostExtend g else deleteCostOpen g - writeArray matrix (ixS, ixT, Insert) $ if predT + insCost == maxScore then insertCostExtend g else insertCostOpen g - writeArray matrix (ixS, ixT, Match) maxScore + | otherwise -> fillInnerMatrix g matrix False sub ixS ixT pure matrix (lowerS, upperS) = bounds s @@ -230,7 +249,7 @@ -- Conditions of traceback are described below -- {-# INLINE cond #-} - cond (LocalAlignment subC gap) = Conditions localStop (defDiag subC) (vert gap) (horiz gap) + cond (LocalAlignment _ gap) = Conditions localStop (move gap) -- Start from bottom right corner -- @@ -280,23 +299,7 @@ writeArray matrix (ixS, ixT, Match) 0 writeArray matrix (ixS, ixT, Insert) $ insertCostOpen g writeArray matrix (ixS, ixT, Delete) $ deleteCostOpen g - | otherwise -> do - predDiag <- matrix `readArray` (pred ixS, pred ixT, Match) - predS <- matrix `readArray` (pred ixS, ixT, Match) - predT <- matrix `readArray` ( ixS, pred ixT, Match) - - delCost <- matrix `readArray` (pred ixS, ixT, Delete) - insCost <- matrix `readArray` ( ixS, pred ixT, Insert) - - let maxScore = maximum [ predDiag + sub ixS ixT - , predS + delCost - , predT + insCost - , 0 - ] - - writeArray matrix (ixS, ixT, Delete) $ if predS + delCost == maxScore then deleteCostExtend g else deleteCostOpen g - writeArray matrix (ixS, ixT, Insert) $ if predT + insCost == maxScore then insertCostExtend g else insertCostOpen g - writeArray matrix (ixS, ixT, Match) maxScore + | otherwise -> fillInnerMatrix g matrix True sub ixS ixT pure matrix (lowerS, upperS) = bounds s @@ -317,7 +320,7 @@ -- Conditions of traceback are described below -- {-# INLINE cond #-} - cond (SemiglobalAlignment subC gap) = Conditions defStop (defDiag subC) (vert gap) (horiz gap) + cond (SemiglobalAlignment _ gap) = Conditions defStop (move gap) -- Start from bottom right corner -- @@ -357,31 +360,11 @@ matrix <- newArray ((lowerS, lowerT, Insert), (nilS, nilT, Match)) 0 :: ST s (STUArray s (Index m, Index m', EditOp) Int) forM_ [lowerS .. nilS] $ \ixS -> forM_ [lowerT .. nilT] $ \ixT -> - - -- Next cell = max (d_i-1,j + gap, d_i,j-1 + gap, d_i-1,j-1 + s(i,j)) - -- Matrices with gap costs are also filled as follows: - -- gepMatrix[i, j] <- gapExtend if one of strings has gap at this position else gapOpen - -- if | ixS == lowerS || ixT == lowerT -> do writeArray matrix (ixS, ixT, Match) 0 writeArray matrix (ixS, ixT, Insert) $ insertCostOpen g writeArray matrix (ixS, ixT, Delete) $ deleteCostOpen g - | otherwise -> do - predDiag <- matrix `readArray` (pred ixS, pred ixT, Match) - predS <- matrix `readArray` (pred ixS, ixT, Match) - predT <- matrix `readArray` ( ixS, pred ixT, Match) - - delCost <- matrix `readArray` (pred ixS, ixT, Delete) - insCost <- matrix `readArray` ( ixS, pred ixT, Insert) - - let maxScore = maximum [ predDiag + sub ixS ixT - , predS + delCost - , predT + insCost - ] - - writeArray matrix (ixS, ixT, Delete) $ if predS + delCost == maxScore then deleteCostExtend g else deleteCostOpen g - writeArray matrix (ixS, ixT, Insert) $ if predT + insCost == maxScore then insertCostExtend g else insertCostOpen g - writeArray matrix (ixS, ixT, Match) maxScore + | otherwise -> fillInnerMatrix g matrix False sub ixS ixT pure matrix (lowerS, upperS) = bounds s @@ -391,3 +374,33 @@ sub :: Index m -> Index m' -> Int sub = substitute subC s t + +{-# INLINE fillInnerMatrix #-} +fillInnerMatrix + :: (IsGap g, Ix ix, Ix ix', Enum ix, Enum ix') + => g + -> STUArray s (ix, ix', EditOp) Int + -> Bool -- ^ Is this local alignment? + -> (ix -> ix' -> Int) -- ^ Substitution function + -> ix + -> ix' + -> ST s () +fillInnerMatrix g matrix isLocal sub ixS ixT = do + predDiag <- matrix `readArray` (pred ixS, pred ixT, Match) + predS <- matrix `readArray` (pred ixS, ixT, Match) + predT <- matrix `readArray` ( ixS, pred ixT, Match) + + delCost <- matrix `readArray` (pred ixS, ixT, Delete) + insCost <- matrix `readArray` ( ixS, pred ixT, Insert) + + let + updInsCost = max (insCost + insertCostExtend g) (predT + insertCostOpen g) + updDelCost = max (delCost + deleteCostExtend g) (predS + deleteCostOpen g) + repCost = predDiag + sub ixS ixT + maxCost = max repCost $ max updInsCost updDelCost + + finalCost = if isLocal then max 0 maxCost else maxCost + + writeArray matrix (ixS, ixT, Delete) updDelCost + writeArray matrix (ixS, ixT, Insert) updInsCost + writeArray matrix (ixS, ixT, Match) finalCost diff --git a/src/Bio/Chain/Alignment/Type.hs b/src/Bio/Chain/Alignment/Type.hs --- a/src/Bio/Chain/Alignment/Type.hs +++ b/src/Bio/Chain/Alignment/Type.hs @@ -111,16 +111,24 @@ -- type Matrix m m' = UArray (Index m, Index m', EditOp) Int --- | Traceback condition type +-- | Traceback stop condition type -- -type Condition m m' = Matrix m m' -> m -> m' -> Index m -> Index m' -> Bool +type Stop m m' = Matrix m m' -> m -> m' -> Index m -> Index m' -> Bool +-- | Traceback next move generator type +-- +type Move m m' + = Matrix m m' + -> m -> m' + -> Index m -> Index m' + -> EditOp -- ^ Current matrix, used in affine alignment + -> (EditOp, Index m, Index m', Operation (Index m) (Index m')) + -- ^ Next matrix, next indices, new operation + -- | A set of traceback conditions -- -data Conditions m m' = Conditions { isStop :: Condition m m' -- ^ Should we stop? - , isDiag :: Condition m m' -- ^ Should we go daigonally? - , isVert :: Condition m m' -- ^ Should we go vertically? - , isHoriz :: Condition m m' -- ^ Should we go horizontally? +data Conditions m m' = Conditions { isStop :: Stop m m' -- ^ Should we stop? + , doMove :: Move m m' -- ^ Where to go next? } -- | Sequence Alignment result diff --git a/test/HandcraftedSpec.hs b/test/HandcraftedSpec.hs --- a/test/HandcraftedSpec.hs +++ b/test/HandcraftedSpec.hs @@ -55,3 +55,17 @@ it "first sequence" $ a4' `shouldBe` a4Ans it "second sequence" $ b4' `shouldBe` b4Ans it "score" $ score4 `shouldBe` 42 + + describe "SemiglobalAlignment. AffineGap. Validate traceback" $ do + let testAlign = align (SemiglobalAlignment blosum62 (AffineGap (-11) (-1))) + let s1 = "SSSLNRTTT" + s2 = "SSSTTT" + s1ans = "SSSLNRTTT" + s2ans = "SSS---TTT" + + res = testAlign s1 s2 + (s1', s2') = viewAlignment res + + it "first sequence" $ s1' `shouldBe` s1ans + it "second sequence" $ s2' `shouldBe` s2ans + it "score" $ score res `shouldBe` 14 diff --git a/test/JuliaSpec.hs b/test/JuliaSpec.hs --- a/test/JuliaSpec.hs +++ b/test/JuliaSpec.hs @@ -78,9 +78,9 @@ res2 = testAlign "AGGTACC" "CAATG" -- Честно посчитано руками в тетради - it "score" $ score res1 `shouldBe` (-4) + it "score" $ score res1 `shouldBe` (-2) it "score" $ score res2 `shouldBe` (-8) - it "alignment" $ viewAlignment res1 `shouldBe` ("-AGGT", "CAA-T") + it "alignment" $ viewAlignment res1 `shouldBe` ("--AGGT", "CAA--T") it "alignment" $ viewAlignment res2 `shouldBe` ("--AGGTACC","CAATG----") describe "GlobalAlignment. Sanity test" $ do