cobot (empty) → 0.1.1.0
raw patch · 32 files changed
+3280/−0 lines, 32 filesdep +arraydep +basedep +bytestringsetup-changed
Dependencies added: array, base, bytestring, clock, cobot, containers, criterion, deepseq, hspec, lens, linear, megaparsec, mtl, parallel, random, split, template-haskell, text
Files
- ChangeLog.md +10/−0
- LICENSE +30/−0
- README.md +4/−0
- Setup.hs +2/−0
- bench/Main.hs +75/−0
- cobot.cabal +130/−0
- src/Bio/Chain.hs +84/−0
- src/Bio/Chain/Alignment.hs +295/−0
- src/Bio/Chain/Alignment/Algorithms.hs +393/−0
- src/Bio/Chain/Alignment/Scoring.hs +109/−0
- src/Bio/Chain/Alignment/Scoring/Loader.hs +30/−0
- src/Bio/Chain/Alignment/Scoring/TH.hs +56/−0
- src/Bio/Chain/Alignment/Type.hs +164/−0
- src/Bio/Molecule.hs +70/−0
- src/Bio/NucleicAcid/Chain.hs +24/−0
- src/Bio/NucleicAcid/Nucleotide.hs +8/−0
- src/Bio/NucleicAcid/Nucleotide/Instances.hs +60/−0
- src/Bio/NucleicAcid/Nucleotide/Type.hs +94/−0
- src/Bio/Protein/Algebra.hs +203/−0
- src/Bio/Protein/AminoAcid.hs +9/−0
- src/Bio/Protein/AminoAcid/Instances.hs +326/−0
- src/Bio/Protein/AminoAcid/Type.hs +299/−0
- src/Bio/Protein/Chain.hs +25/−0
- src/Bio/Protein/Chain/Builder.hs +141/−0
- src/Bio/Protein/Metric.hs +28/−0
- src/Bio/Utils/Geometry.hs +108/−0
- src/Bio/Utils/IUPAC.hs +39/−0
- src/Bio/Utils/Matrix.hs +20/−0
- src/Bio/Utils/Monomer.hs +30/−0
- test/HandcraftedSpec.hs +57/−0
- test/JuliaSpec.hs +270/−0
- test/Spec.hs +87/−0
+ ChangeLog.md view
@@ -0,0 +1,10 @@+# Changelog for cobot++## Unreleased changes++## [0.1.1.0] - 2019-06-17+### Added+- Typeclass `IsGap`.+- Now you can align two sequences using different gaps for each one of them.+### Changed+- Removed `affine` function from `SequenceAlignment` typeclass.
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright Pavel Yakovlev (c) 2018++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 Pavel Yakovlev 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.
+ README.md view
@@ -0,0 +1,4 @@+# cobot++[](https://travis-ci.org/less-wrong/cobot)+[]()
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ bench/Main.hs view
@@ -0,0 +1,75 @@+module Main where++import Bio.Chain (Chain, fromList)+import Bio.Chain.Alignment (AffineGap (..),+ GlobalAlignment (..),+ LocalAlignment (..),+ SemiglobalAlignment (..),+ SimpleGap, align)+import Bio.Chain.Alignment.Scoring (nuc44)+import Control.DeepSeq (NFData (..), deepseq)+import Control.Monad (replicateM)+import Control.Monad.State (State, evalState, state)+import Control.Parallel.Strategies (dot, parListChunk, rdeepseq, rpar,+ withStrategy)+import Criterion (bench, bgroup, env, nfIO)+import Criterion.Main (defaultMain)+import GHC.Conc (numCapabilities)+import System.Clock (Clock (Monotonic), diffTimeSpec,+ getTime, nsec, sec)+import System.Random (RandomGen, getStdGen, randomR)++makeRandomChain :: RandomGen g => Int -> State g String+makeRandomChain 0 = pure ""+makeRandomChain len = do+ c <- ("ATGC" !!) <$> state (randomR (0, 3))+ cs <- makeRandomChain (len - 1)+ pure (c : cs)++makeRandomChainIO :: Int -> IO (Chain Int Char)+makeRandomChainIO len = do+ list <- evalState (makeRandomChain len) <$> getStdGen+ pure (fromList list)++measureTime :: NFData a => String -> a -> IO ()+measureTime label value = do+ t1 <- getTime Monotonic+ t2 <- value `deepseq` getTime Monotonic+ let dt = diffTimeSpec t2 t1+ let timeInSeconds = fromIntegral (sec dt) + fromIntegral (nsec dt) / 1000000000 :: Double+ let padding = " " <> replicate (max 0 (50 - length label)) '.' <> " "+ putStrLn $ label <> padding <> show timeInSeconds <> "s"++parMap' :: NFData b => Int -> (a -> b) -> [a] -> [b]+parMap' chunkSize f = withStrategy (parListChunk chunkSize (rdeepseq `dot` rpar)) . map f++setupEnv :: IO (Chain Int Char, [Chain Int Char], Int)+setupEnv = do+ a <- makeRandomChainIO 600+ bs <- replicateM 20 $ makeRandomChainIO 4500+ let chunkSize = length bs `div` numCapabilities+ pure (a, bs, chunkSize)++main :: IO ()+main = defaultMain [+ env setupEnv $ \ ~(a, bs, chunkSize) -> bgroup "main" [+ bench "Local alignment" $+ let align' = align (LocalAlignment nuc44 (-10 :: SimpleGap)) a+ in nfIO . pure $ parMap' chunkSize align' bs,+ bench "Global alignment" $+ let align' = align (GlobalAlignment nuc44 (-10 :: SimpleGap)) a+ in nfIO . pure $ parMap' chunkSize align' bs,+ bench "Semiglobal alignment" $+ let align' = align (SemiglobalAlignment nuc44 (-10 :: SimpleGap)) a+ in nfIO . pure $ parMap' chunkSize align' bs,+ bench "Local alignment with affine gap" $+ let align' = align (LocalAlignment nuc44 (AffineGap (-10) (-1))) a+ in nfIO . pure $ parMap' chunkSize align' bs,+ bench "Global alignment with affine gap" $+ let align' = align (GlobalAlignment nuc44 (AffineGap (-10) (-1))) a+ in nfIO . pure $ parMap' chunkSize align' bs,+ bench "Semiglobal alignment with affine gap" $+ let align' = align (SemiglobalAlignment nuc44 (AffineGap (-10) (-1))) a+ in nfIO . pure $ parMap' chunkSize align' bs+ ]+ ]
+ cobot.cabal view
@@ -0,0 +1,130 @@+cabal-version: 1.12++-- This file has been generated from package.yaml by hpack version 0.31.2.+--+-- see: https://github.com/sol/hpack+--+-- hash: 5253d8bdc8faf78e15758bbf5ca466c06c743718a111840ff74966c4932315ad++name: cobot+version: 0.1.1.0+synopsis: Computational biology toolkit to collaborate with researchers in constructive protein engineering+description: Please see the README on GitHub at <https://github.com/less-wrong/cobot#readme>+category: Bio+homepage: https://github.com/less-wrong/cobot#readme+bug-reports: https://github.com/less-wrong/cobot/issues+author: Pavel Yakovlev, Bogdan Neterebskii, Alexander Sadovnikov+maintainer: pavel@yakovlev.me+copyright: 2018—2019, Less Wrong Bio+license: BSD3+license-file: LICENSE+build-type: Simple+extra-source-files:+ README.md+ ChangeLog.md++source-repository head+ type: git+ location: https://github.com/less-wrong/cobot++library+ exposed-modules:+ Bio.Chain+ Bio.Chain.Alignment+ Bio.Chain.Alignment.Algorithms+ Bio.Chain.Alignment.Scoring+ Bio.Chain.Alignment.Scoring.Loader+ Bio.Chain.Alignment.Scoring.TH+ Bio.Chain.Alignment.Type+ Bio.Molecule+ Bio.NucleicAcid.Chain+ Bio.NucleicAcid.Nucleotide+ Bio.NucleicAcid.Nucleotide.Instances+ Bio.NucleicAcid.Nucleotide.Type+ Bio.Protein.Algebra+ Bio.Protein.AminoAcid+ Bio.Protein.AminoAcid.Instances+ Bio.Protein.AminoAcid.Type+ Bio.Protein.Chain+ Bio.Protein.Chain.Builder+ Bio.Protein.Metric+ Bio.Utils.Geometry+ Bio.Utils.IUPAC+ Bio.Utils.Matrix+ Bio.Utils.Monomer+ other-modules:+ Paths_cobot+ hs-source-dirs:+ src+ default-extensions: AllowAmbiguousTypes ConstraintKinds DeriveFoldable DeriveFunctor DeriveGeneric DeriveTraversable FlexibleContexts FlexibleInstances GeneralizedNewtypeDeriving MultiWayIf RankNTypes RecordWildCards ScopedTypeVariables TypeApplications TypeFamilies TypeSynonymInstances UndecidableInstances+ build-depends:+ array+ , base >=4.7 && <5+ , bytestring+ , containers+ , deepseq+ , lens+ , linear+ , megaparsec+ , mtl+ , split+ , template-haskell+ , text+ default-language: Haskell2010++test-suite cobot-test+ type: exitcode-stdio-1.0+ main-is: Spec.hs+ other-modules:+ HandcraftedSpec+ JuliaSpec+ Paths_cobot+ hs-source-dirs:+ test+ default-extensions: AllowAmbiguousTypes ConstraintKinds DeriveFoldable DeriveFunctor DeriveGeneric DeriveTraversable FlexibleContexts FlexibleInstances GeneralizedNewtypeDeriving MultiWayIf RankNTypes RecordWildCards ScopedTypeVariables TypeApplications TypeFamilies TypeSynonymInstances UndecidableInstances+ ghc-options: -threaded -rtsopts "-with-rtsopts=-A64m -qb0 -I0 -N -qn4"+ build-depends:+ array+ , base >=4.7 && <5+ , bytestring+ , cobot+ , containers+ , deepseq+ , hspec+ , lens+ , linear+ , megaparsec+ , mtl+ , split+ , template-haskell+ , text+ default-language: Haskell2010++benchmark cobot-bench+ type: exitcode-stdio-1.0+ main-is: Main.hs+ other-modules:+ Paths_cobot+ hs-source-dirs:+ bench+ default-extensions: AllowAmbiguousTypes ConstraintKinds DeriveFoldable DeriveFunctor DeriveGeneric DeriveTraversable FlexibleContexts FlexibleInstances GeneralizedNewtypeDeriving MultiWayIf RankNTypes RecordWildCards ScopedTypeVariables TypeApplications TypeFamilies TypeSynonymInstances UndecidableInstances OverloadedStrings+ ghc-options: -threaded -rtsopts "-with-rtsopts=-A64m -qb0 -I0 -N -qn4"+ build-depends:+ array+ , base >=4.7 && <5+ , bytestring+ , clock+ , cobot+ , containers+ , criterion+ , deepseq+ , lens+ , linear+ , megaparsec+ , mtl+ , parallel+ , random+ , split+ , template-haskell+ , text+ default-language: Haskell2010
+ src/Bio/Chain.hs view
@@ -0,0 +1,84 @@+{-# LANGUAGE TupleSections #-}+module Bio.Chain+ ( ChainLike (..)+ , Chain+ , chain, fromList+ , (!), (//)+ ) where++import Control.Lens+import qualified Data.Array as A ( bounds+ , assocs+ )+import Data.Array ( Array+ , Ix+ , array+ , listArray+ , (!)+ , (//)+ )+import Data.Array.Base (unsafeAt)++type Chain i a = Array i a++-- | Construct new chain from list+--+chain :: Ix i => (i, i) -> [(i, a)] -> Chain i a+chain = array++-- | Construct new int-labeled chain from list+--+fromList :: [a] -> Chain Int a+fromList lst = listArray (0, length lst - 1) lst++-- | Chain-like sequence, by default it is an array or a list+--+class (Ixed m, Enum (Index m)) => ChainLike m where+ bounds :: m -> (Index m, Index m)+ assocs :: m -> [(Index m, IxValue m)]+ modify :: Index m -> (IxValue m -> IxValue m) -> m -> m+ modifyBefore :: Index m -> (IxValue m -> IxValue m) -> m -> m+ modifyAfter :: Index m -> (IxValue m -> IxValue m) -> m -> m++ unsafeRead :: m -> Index m -> IxValue m+ unsafeRead ch i = ch ^?! ix i++instance ChainLike [a] where+ bounds = (0,) . pred . length++ assocs = zip [0..]++ modify _ _ [] = []+ modify 0 f (x:xs) = f x:xs+ modify i f (x:xs) = x:modify (i - 1) f xs++ modifyBefore i f lst = (f <$> take i lst) ++ drop i lst+ modifyAfter i f lst = take (i + 1) lst ++ (f <$> drop (i + 1) lst)++ unsafeRead = (!!)++instance (Ix i, Enum i) => ChainLike (Array i a) where+ bounds = A.bounds++ assocs = A.assocs++ modify i f ar = ar // [(i, f (ar ! i))]++ modifyBefore i f ar = let (mi, _) = bounds ar+ in ar // [(j, f (ar ! j)) | j <- [mi .. pred i]]+ modifyAfter i f ar = let (_, ma) = bounds ar+ in ar // [(j, f (ar ! j)) | j <- [succ i .. ma]]++ {-# INLINE unsafeRead #-}+ unsafeRead = unsafeReadArray++class (Ixed m) => UnsafeReadArray m where+ unsafeReadArray :: m -> Index m -> IxValue m++instance (Ix i, Enum i) => UnsafeReadArray (Array i a) where+ {-# INLINE unsafeReadArray #-}+ unsafeReadArray = (!)++instance {-# OVERLAPPING #-} UnsafeReadArray (Array Int a) where+ {-# INLINE unsafeReadArray #-}+ unsafeReadArray = unsafeAt
+ src/Bio/Chain/Alignment.hs view
@@ -0,0 +1,295 @@+module Bio.Chain.Alignment+ ( AlignmentResult (..), SimpleGap, SimpleGap2, AffineGap (..), AffineGap2, Operation (..)+ , EditDistance (..)+ , GlobalAlignment (..), LocalAlignment (..), SemiglobalAlignment (..)+ , IsGap (..)+ , align+ , viewAlignment+ , prettyAlignmment+ , similarityGen+ , differenceGen+ , similarity+ , difference+ ) where++import Control.Lens (Index, IxValue, Ixed (..), to,+ (^?!))+import Data.Array.Unboxed ((!))+import Data.List (intercalate)+import Data.List.Split (chunksOf)++import Bio.Chain hiding ((!))+import Bio.Chain.Alignment.Algorithms+import Bio.Chain.Alignment.Type+import Bio.Utils.Geometry (R)+import Bio.Utils.Monomer (Symbol (..))++-- | Align chains using specifed algorithm+--+{-# SPECIALISE align :: LocalAlignment SimpleGap Char Char -> Chain Int Char -> Chain Int Char -> AlignmentResult (Chain Int Char) (Chain Int Char) #-}+{-# SPECIALISE align :: LocalAlignment AffineGap Char Char -> Chain Int Char -> Chain Int Char -> AlignmentResult (Chain Int Char) (Chain Int Char) #-}+{-# SPECIALISE align :: SemiglobalAlignment SimpleGap Char Char -> Chain Int Char -> Chain Int Char -> AlignmentResult (Chain Int Char) (Chain Int Char) #-}+{-# SPECIALISE align :: SemiglobalAlignment AffineGap Char Char -> Chain Int Char -> Chain Int Char -> AlignmentResult (Chain Int Char) (Chain Int Char) #-}+{-# SPECIALISE align :: GlobalAlignment SimpleGap Char Char -> Chain Int Char -> Chain Int Char -> AlignmentResult (Chain Int Char) (Chain Int Char) #-}+{-# SPECIALISE align :: GlobalAlignment AffineGap Char Char -> Chain Int Char -> Chain Int Char -> AlignmentResult (Chain Int Char) (Chain Int Char) #-}+align :: forall algo m m'.(SequenceAlignment algo, Alignable m, Alignable m') => algo (IxValue m) (IxValue m') -> m -> m' -> AlignmentResult m m'+align algo s t = AlignmentResult alignmentScore alignmentResult s t+ where+ -- Bounds of chains specify bounds of alignment matrix+ (lowerS, upperS) = bounds s+ (lowerT, upperT) = bounds t+ -- Fill the matrix+ mat :: Matrix m m'+ mat = scoreMatrix algo s t+ -- Result coordinates+ coords :: (Index m, Index m')+ coords = traceStart algo mat s t+ -- Score of alignment+ alignmentScore :: Int+ alignmentScore = let (x, y) = coords in mat ! (x, y, Match)++ -- Resulting alignment should contain additional deletions/insertions in case of semiglobal+ -- alignment+ alignmentResult :: [Operation (Index m) (Index m')]+ alignmentResult+ | semi algo = preResult ++ suffix+ | otherwise = preResult+ where+ preResult = uncurry (traceback algo mat s t) coords+ -- Last index of FIRST chain affected by some operation in preResult or (lowerS - 1).+ lastI = last . (pred lowerS :) . map getI $ filter (not . isInsert) preResult+ -- Last index of SECOND chain affected by some operation in preResult or (lowerS - 1).+ lastJ = last . (pred lowerT :) . map getJ $ filter (not . isDelete) preResult+ -- Deletions and insertions of symbols after last operation in preResult+ suffix = case last (MATCH (pred lowerS) (pred lowerT) : preResult) of+ MATCH i j -> map DELETE [succ i .. upperS] ++ map INSERT [succ j .. upperT]+ INSERT _ -> map DELETE [succ lastI .. upperS]+ DELETE _ -> map INSERT [succ lastJ .. upperT]++-- | Traceback function.+--+-- Builds traceback for alignment algorithm @algo@ in matrix @mat@, that is+-- result of alignment of sequences @s@ and @t@. Traceback will start from+-- position with coordinates (@i@, @j@) in matrix.+--+-- Traceback is represented as list of 'Operation's.+--+traceback :: (SequenceAlignment algo, Alignable m, Alignable m')+ => algo (IxValue m) (IxValue m')+ -> Matrix m m'+ -> m+ -> m'+ -> Index m+ -> Index m'+ -> [Operation (Index m) (Index m')]+traceback algo mat s t i' j' = helper i' j' []+ 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"++---------------------------------------------------------------------------------------------------------+ --+ -- 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")+ --+---------------------------------------------------------------------------------------------------------++-- | Calculate similarity and difference between two sequences, aligning them first using given algorithm.+--+similarityGen :: forall algo m m'.(SequenceAlignment algo, Alignable m, Alignable m')+ => algo (IxValue m) (IxValue m')+ -> (IxValue m -> IxValue m' -> Bool)+ -> m+ -> m'+ -> R+similarityGen algo genericEq s t = fromIntegral hamming / fromIntegral len+ where+ operations = alignment (align algo s t)+ len = length operations+ hamming = sum $ toScores <$> operations++ toScores :: Operation (Index m) (Index m') -> Int+ toScores (MATCH i j) = if (s ^?! ix i) `genericEq` (t ^?! ix j) then 1 else 0+ toScores _ = 0++similarity :: 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')+ -> m+ -> m'+ -> R+similarity algo = similarityGen algo (==)+++differenceGen :: forall algo m m'.(SequenceAlignment algo, Alignable m, Alignable m')+ => algo (IxValue m) (IxValue m')+ -> (IxValue m -> IxValue m' -> Bool)+ -> m+ -> m'+ -> R+differenceGen algo genericEq s t = 1.0 - similarityGen algo genericEq s t+++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')+ -> m+ -> m'+ -> R+difference algo = differenceGen algo (==)++-- | 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)+ where+ (s, t) = (sequence1 ar, sequence2 ar)++ toChars :: Operation (Index m) (Index m') -> (Char, Char)+ toChars (MATCH i j) = (symbol (s ^?! ix i), symbol (t ^?! ix j))+ toChars (DELETE i) = (symbol (s ^?! ix i), '-')+ toChars (INSERT j) = ('-', symbol (t ^?! ix j))++-- | Format alignment result as pretty columns of symbols.+--+-- Example with width equal to 20:+--+-- @+-- 0 -------------------- 0+--+-- 0 TTTTTTTTTTTTTTTTTTTT 19+--+-- 0 --GCCTGAATGGTGTGGTGT 17+-- || |||||| |||| |||+-- 20 TTGC-TGAATG-TGTG-TGT 36+--+-- 18 TCGGCGGAGGGACCCAGCTA 37+-- || |||||||||||||||+-- 37 -CG-CGGAGGGACCCAGCT- 53+--+-- 38 AAAAAAAAAA 47+--+-- 53 ---------- 53+-- @+prettyAlignmment+ :: forall m m'+ . (Alignable m, Alignable m', Symbol (IxValue m), Symbol (IxValue m'))+ => AlignmentResult m m' -- ^ Result of alignment to format+ -> Int -- ^ Desired width of one alignment row+ -> String+prettyAlignmment ar width =+ -- Due to construction 'resultRows' first element will be empty string, and we don't need it.+ intercalate "\n" $ tail resultRows+ where+ (s, t) = (sequence1 ar, sequence2 ar)+ rows = chunksOf width $ alignment ar++ chainLength :: forall c. ChainLike c => c -> Int+ chainLength ch = let (a, b) = bounds ch in fromEnum b - fromEnum a + 1++ -- Determine how many characters to leave for position numbers+ numWidth = length $ show $ max (chainLength s) (chainLength t)++ padLeft :: String -> String+ padLeft x = replicate (numWidth - length x) ' ' <> x++ -- Build one column of nice alignment like+ -- T+ -- |+ -- T+ toCharTriple :: Operation (Index m) (Index m') -> (Char, Char, Char)+ toCharTriple (MATCH i j) = (left, if left == right then '|' else ' ', right)+ where+ left = s ^?! ix i . to symbol+ right = t ^?! ix j . to symbol+ toCharTriple (DELETE i) = (s ^?! ix i . to symbol, ' ', '-')+ toCharTriple (INSERT j) = ('-', ' ', t ^?! ix j . to symbol)++ -- Format one chunk of alignment, adding indices to start and end of strings+ -- (prevI, prevJ) must be 1-based indices of last printed characters in each string.+ formatRow :: (Int, Int) -> [Operation (Index m) (Index m')] -> ((Int, Int), [String])+ formatRow (prevI, prevJ) row = ((lastI, lastJ), [resLine1, resLine2, resLine3])+ where+ (line1, line2, line3) = unzip3 $ map toCharTriple row++ -- | Folding function to count lengths of both strings in alignment row+ countChars :: (Int, Int) -> Operation (Index m) (Index m') -> (Int, Int)+ countChars (li, lj) (MATCH _ _) = (li + 1, lj + 1)+ countChars (li, lj) (DELETE _) = (li + 1, lj)+ countChars (li, lj) (INSERT _) = (li, lj + 1)++ (lengthI, lengthJ) = foldl countChars (0, 0) row++ -- Indices of first printed non-gap characters in the current row.+ -- If the row contains any non-gap characters, this is equal to index+ -- of last printed character + 1+ (firstI, firstJ) =+ ( if lengthI > 0 then prevI + 1 else prevI+ , if lengthJ > 0 then prevJ + 1 else prevJ+ )+ (lastI, lastJ) = (prevI + lengthI, prevJ + lengthJ)++ -- It's easier to do everything in 1-based indices and convert before showing+ toZeroBased :: Int -> Int+ toZeroBased 0 = 0+ toZeroBased i = i - 1++ resLine1 = padLeft (show $ toZeroBased firstI) <> " " <> line1 <> " " <> padLeft (show $ toZeroBased lastI)+ resLine2 = padLeft "" <> " " <> line2+ resLine3 = padLeft (show $ toZeroBased firstJ) <> " " <> line3 <> " " <> padLeft (show $ toZeroBased lastJ)++ -- Go through all chunks of operations, accummulating current offsets in both strings+ (_, resultRows) =+ foldl+ (\(off, res) ops -> let (newOff, newRes) = formatRow off ops in (newOff, res <> [""] <> newRes))+ ((0, 0), [])+ rows
+ src/Bio/Chain/Alignment/Algorithms.hs view
@@ -0,0 +1,393 @@+{-# LANGUAGE InstanceSigs #-}+{-# LANGUAGE TupleSections #-}+module Bio.Chain.Alignment.Algorithms where++import Control.Lens (Index, IxValue)+import qualified Data.Array.Unboxed as A (bounds, range)+import Data.List (maximumBy)++import Bio.Chain hiding ((!))+import Bio.Chain.Alignment.Type+import Control.Monad (forM_)+import Data.Ord (comparing)++import Control.Monad.ST (ST)+import Data.Array.Base (readArray, writeArray)+import Data.Array.ST (MArray (..), STUArray, newArray,+ runSTUArray)+import Data.Array.Unboxed (Ix (..), UArray, (!))+++-- | 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+data SemiglobalAlignment a e1 e2 = SemiglobalAlignment (Scoring e1 e2) a++-- Common functions++-- | Lift simple substitution function to a ChainLike collection+--+{-# SPECIALISE substitute :: (Char -> Char -> Int) -> Chain Int Char -> Chain Int Char -> Int -> Int -> Int #-}+{-# INLINE substitute #-}+substitute :: (Alignable m, Alignable m') => (IxValue m -> IxValue m' -> Int) -> m -> m' -> Index m -> Index m' -> Int+substitute f s t i j = f (s `unsafeRead` (pred i)) (t `unsafeRead` (pred j))++-- | Simple substitution function for edit distance+--+substituteED :: EditDistance e1 e2 -> (e1 -> e2 -> Int)+substituteED (EditDistance genericEq) x y = if x `genericEq` y then 1 else 0++-- | Default traceback stop condition.+--+{-# SPECIALISE defStop :: Matrix (Chain Int Char) (Chain Int Char) -> Chain Int Char -> Chain Int Char -> Int -> Int -> Bool #-}+{-# INLINE defStop #-}+defStop :: (Alignable m, Alignable m') => Matrix m m' -> m -> m' -> Index m -> Index m' -> Bool+defStop _ s t i j = let (lowerS, _) = bounds s+ (lowerT, _) = bounds t+ in i == lowerS && j == lowerT++-- | Traceback stop condition for the local alignment.+--+{-# SPECIALISE localStop :: Matrix (Chain Int Char) (Chain Int Char) -> Chain Int Char -> Chain Int Char -> Int -> Int -> Bool #-}+{-# INLINE localStop #-}+localStop :: (Alignable m, Alignable m') => Matrix m m' -> m -> m' -> Index m -> Index m' -> Bool+localStop m' s t i j = let (lowerS, _) = bounds s+ (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++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++ (lowerT, _) = bounds t+ (lowerS, _) = bounds s++-- | Default condition of moving diagonally in traceback.+--+{-# SPECIALISE defDiag :: (Char -> Char -> Int) -> Matrix (Chain Int Char) (Chain Int Char) -> Chain Int Char -> Chain Int Char -> Int -> Int -> Bool #-}+{-# INLINE defDiag #-}+defDiag :: (Alignable m, Alignable m') => (IxValue m -> IxValue m' -> Int) -> Matrix m m' -> m -> m' -> Index m -> Index m' -> Bool+defDiag sub' m s t i j = let sub = substitute sub' s t+ in m ! (pred i, pred j, Match) + sub i j == m ! (i, j, Match)++-- | Default start condition for traceback.+--+{-# SPECIALISE defStart :: Matrix (Chain Int Char) (Chain Int Char) -> Chain Int Char -> Chain Int Char -> (Int, Int) #-}+{-# INLINE defStart #-}+defStart :: (Alignable m, Alignable m') => Matrix m m' -> m -> m' -> (Index m, Index m')+defStart m _ _ = let ((_, _, _), (upperS, upperT, _)) = A.bounds m in (upperS, upperT)++-- | Default start condition for traceback in local alignment.+--+{-# SPECIALISE localStart :: Matrix (Chain Int Char) (Chain Int Char) -> Chain Int Char -> Chain Int Char -> (Int, Int) #-}+{-# INLINE localStart #-}+localStart :: (Alignable m, Alignable m') => Matrix m m' -> m -> m' -> (Index m, Index m')+localStart m _ _ = let ((lowerS, lowerT, _), (upperS, upperT, _)) = A.bounds m+ range' = A.range ((lowerS, lowerT, Match), (upperS, upperT, Match))+ in (\(a, b, _) -> (a, b)) $ maximumBy (comparing (m !)) range'++-- | Default start condition for traceback in semiglobal alignment.+--+{-# SPECIALISE semiStart :: Matrix (Chain Int Char) (Chain Int Char) -> Chain Int Char -> Chain Int Char -> (Int, Int) #-}+{-# INLINE semiStart #-}+semiStart :: (Alignable m, Alignable m') => Matrix m m' -> m -> m' -> (Index m, Index m')+semiStart m _ _ = let ((lowerS, lowerT, _), (upperS, upperT, _)) = A.bounds m+ lastCol = (, upperT, Match) <$> [lowerS .. upperS]+ lastRow = (upperS, , Match) <$> [lowerT .. upperT]+ in (\(a, b, _) -> (a, b)) $ maximumBy (comparing (m !)) $ lastCol ++ lastRow++-- Alignment algorithm instances++-------------------+ --+ -- About affine gaps+ --+ -- There are three matrices used in all the algorithms below:+ -- 1) One stores the resulting scores for each prefix pair;+ -- 2) One stores insertion costs in the first sequence for each prefix pair;+ -- 3) One stores insertion costs in the second sequence for each prefix pair.+ --+ -- Matrices 2 and 3 are used in affine penalty calculation:+ -- Let `gapOpen` and `gapExtend` be affine gap penalties (for opening and extending a gap correspondingly).+ -- M2[i, j] = gapOpen if neither of sequences has insertion or deletion at prefix (i, j);+ -- M2[i, j] = gapExtend if sequence1 has insertion or, equivalently, sequence2 has deletion at prefix (i, j);+ --+ -- gap penalty in the first sequence for the prefix (i, j) is gapOpen + M2[i, j]+ -- So, if there are no gaps in the sequence before, the penalty will be `gapOpen`.+ -- Otherwise it will be `gapExtend`.+ --+ -- So, M2 holds insertion costs for the first sequence, and M3 holds insertion costs for the second sequence.+ --+ -- The resulting score is the same as in plain gap penalty:+ -- the biggest one between substitution, insertion and deletion scores.+ --+-------------------++instance IsGap g => SequenceAlignment (GlobalAlignment g) where++ -- Conditions of traceback are described below+ --+ {-# INLINE cond #-}+ cond (GlobalAlignment subC gap) = Conditions defStop (defDiag subC) (vert gap) (horiz gap)++ -- Start from bottom right corner+ --+ {-# INLINE traceStart #-}+ traceStart = const defStart++ scoreMatrix :: forall m m' . (Alignable m, Alignable m')+ => GlobalAlignment g (IxValue m) (IxValue m')+ -> m+ -> m'+ -> Matrix m m'+ scoreMatrix (GlobalAlignment subC g) s t | isAffine g = uMatrixAffine+ | otherwise = uMatrixSimple+ where+ uMatrixSimple :: UArray (Index m, Index m', EditOp) Int+ uMatrixSimple = runSTUArray $ do+ matrix <- newArray ((lowerS, lowerT, Match), (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))+ --+ if | ixS == lowerS -> writeArray matrix (ixS, ixT, Match) $ (insertCostOpen g) * index (lowerT, nilT) ixT+ | ixT == lowerT -> writeArray matrix (ixS, ixT, Match) $ (deleteCostOpen g) * index (lowerS, nilS) ixS+ | otherwise -> do+ predDiag <- matrix `readArray` (pred ixS, pred ixT, Match)+ predS <- matrix `readArray` (pred ixS, ixT, Match)+ predT <- matrix `readArray` ( ixS, pred ixT, Match)+ writeArray matrix (ixS, ixT, Match) $ maximum [ predDiag + sub ixS ixT+ , predS + deleteCostOpen g+ , predT + insertCostOpen g+ ]+ pure matrix++ uMatrixAffine :: UArray (Index m, Index m', EditOp) Int+ uMatrixAffine = runSTUArray $ do+ 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+ | ixS == lowerS -> do+ writeArray matrix (ixS, ixT, Match) $ insertCostOpen g + (insertCostExtend g) * pred (index (lowerT, nilT) ixT)+ writeArray matrix (ixS, ixT, Insert) $ insertCostExtend g+ writeArray matrix (ixS, ixT, Delete) $ deleteCostOpen g+ | ixT == lowerT -> do+ 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+ pure matrix++ (lowerS, upperS) = bounds s+ (lowerT, upperT) = bounds t+ nilS = succ upperS+ nilT = succ upperT++ sub :: Index m -> Index m' -> Int+ sub = substitute subC s t++instance IsGap g => SequenceAlignment (LocalAlignment g) where++ -- Conditions of traceback are described below+ --+ {-# INLINE cond #-}+ cond (LocalAlignment subC gap) = Conditions localStop (defDiag subC) (vert gap) (horiz gap)++ -- Start from bottom right corner+ --+ {-# INLINE traceStart #-}+ traceStart = const localStart++ scoreMatrix :: forall m m' . (Alignable m, Alignable m')+ => LocalAlignment g (IxValue m) (IxValue m')+ -> m+ -> m'+ -> Matrix m m'+ scoreMatrix (LocalAlignment subC g) s t | isAffine g = uMatrixAffine+ | otherwise = uMatrixSimple+ where+ uMatrixSimple :: UArray (Index m, Index m', EditOp) Int+ uMatrixSimple = runSTUArray $ do+ matrix <- newArray ((lowerS, lowerT, Match), (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))+ --+ if | ixS == lowerS -> writeArray matrix (ixS, ixT, Match) 0+ | ixT == lowerT -> writeArray matrix (ixS, ixT, Match) 0+ | otherwise -> do+ predDiag <- matrix `readArray` (pred ixS, pred ixT, Match)+ predS <- matrix `readArray` (pred ixS, ixT, Match)+ predT <- matrix `readArray` ( ixS, pred ixT, Match)+ writeArray matrix (ixS, ixT, Match) $ maximum [ predDiag + sub ixS ixT+ , predS + deleteCostOpen g+ , predT + insertCostOpen g+ , 0+ ]+ pure matrix++ uMatrixAffine :: UArray (Index m, Index m', EditOp) Int+ uMatrixAffine = runSTUArray $ do+ 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+ , 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+ pure matrix++ (lowerS, upperS) = bounds s+ (lowerT, upperT) = bounds t+ nilS = succ upperS+ nilT = succ upperT++ sub :: Index m -> Index m' -> Int+ sub = substitute subC s t++instance IsGap g => SequenceAlignment (SemiglobalAlignment g) where++ -- The alignment is semiglobal, so we have to perform some additional operations+ --+ {-# INLINE semi #-}+ semi = const True++ -- Conditions of traceback are described below+ --+ {-# INLINE cond #-}+ cond (SemiglobalAlignment subC gap) = Conditions defStop (defDiag subC) (vert gap) (horiz gap)++ -- Start from bottom right corner+ --+ {-# INLINE traceStart #-}+ traceStart = const semiStart++ scoreMatrix :: forall m m' . (Alignable m, Alignable m')+ => SemiglobalAlignment g (IxValue m) (IxValue m')+ -> m+ -> m'+ -> Matrix m m'+ scoreMatrix (SemiglobalAlignment subC g) s t | isAffine g = uMatrixAffine+ | otherwise = uMatrixSimple+ where+ uMatrixSimple :: UArray (Index m, Index m', EditOp) Int+ uMatrixSimple = runSTUArray $ do+ matrix <- newArray ((lowerS, lowerT, Match), (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))+ --+ if | ixS == lowerS -> writeArray matrix (ixS, ixT, Match) 0+ | ixT == lowerT -> writeArray matrix (ixS, ixT, Match) 0+ | otherwise -> do+ predDiag <- matrix `readArray` (pred ixS, pred ixT, Match)+ predS <- matrix `readArray` (pred ixS, ixT, Match)+ predT <- matrix `readArray` ( ixS, pred ixT, Match)+ writeArray matrix (ixS, ixT, Match) $ maximum [ predDiag + sub ixS ixT+ , predS + deleteCostOpen g+ , predT + insertCostOpen g+ ]+ pure matrix++ uMatrixAffine :: UArray (Index m, Index m', EditOp) Int+ uMatrixAffine = runSTUArray $ do+ 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+ pure matrix++ (lowerS, upperS) = bounds s+ (lowerT, upperT) = bounds t+ nilS = succ upperS+ nilT = succ upperT++ sub :: Index m -> Index m' -> Int+ sub = substitute subC s t
+ src/Bio/Chain/Alignment/Scoring.hs view
@@ -0,0 +1,109 @@+{-# LANGUAGE QuasiQuotes #-}+{-# LANGUAGE TemplateHaskell #-}+module Bio.Chain.Alignment.Scoring+ ( module L+ , matrix+ , BLOSUM62 (..), PAM250 (..), NUC44 (..)+ , blosum62, pam250, nuc44+ ) where++import Bio.Chain.Alignment.Scoring.Loader as L+import Bio.Chain.Alignment.Scoring.TH ( matrix )++[matrix|BLOSUM62+# Matrix made by matblas from blosum62.iij+# * column uses minimum score+# BLOSUM Clustered Scoring Matrix in 1/2 Bit Units+# Blocks Database = /data/blocks_5.0/blocks.dat+# Cluster Percentage: >= 62+# Entropy = 0.6979, Expected = -0.5209+ A R N D C Q E G H I L K M F P S T W Y V B Z X *+A 4 -1 -2 -2 0 -1 -1 0 -2 -1 -1 -1 -1 -2 -1 1 0 -3 -2 0 -2 -1 0 -4+R -1 5 0 -2 -3 1 0 -2 0 -3 -2 2 -1 -3 -2 -1 -1 -3 -2 -3 -1 0 -1 -4+N -2 0 6 1 -3 0 0 0 1 -3 -3 0 -2 -3 -2 1 0 -4 -2 -3 3 0 -1 -4+D -2 -2 1 6 -3 0 2 -1 -1 -3 -4 -1 -3 -3 -1 0 -1 -4 -3 -3 4 1 -1 -4+C 0 -3 -3 -3 9 -3 -4 -3 -3 -1 -1 -3 -1 -2 -3 -1 -1 -2 -2 -1 -3 -3 -2 -4+Q -1 1 0 0 -3 5 2 -2 0 -3 -2 1 0 -3 -1 0 -1 -2 -1 -2 0 3 -1 -4+E -1 0 0 2 -4 2 5 -2 0 -3 -3 1 -2 -3 -1 0 -1 -3 -2 -2 1 4 -1 -4+G 0 -2 0 -1 -3 -2 -2 6 -2 -4 -4 -2 -3 -3 -2 0 -2 -2 -3 -3 -1 -2 -1 -4+H -2 0 1 -1 -3 0 0 -2 8 -3 -3 -1 -2 -1 -2 -1 -2 -2 2 -3 0 0 -1 -4+I -1 -3 -3 -3 -1 -3 -3 -4 -3 4 2 -3 1 0 -3 -2 -1 -3 -1 3 -3 -3 -1 -4+L -1 -2 -3 -4 -1 -2 -3 -4 -3 2 4 -2 2 0 -3 -2 -1 -2 -1 1 -4 -3 -1 -4+K -1 2 0 -1 -3 1 1 -2 -1 -3 -2 5 -1 -3 -1 0 -1 -3 -2 -2 0 1 -1 -4+M -1 -1 -2 -3 -1 0 -2 -3 -2 1 2 -1 5 0 -2 -1 -1 -1 -1 1 -3 -1 -1 -4+F -2 -3 -3 -3 -2 -3 -3 -3 -1 0 0 -3 0 6 -4 -2 -2 1 3 -1 -3 -3 -1 -4+P -1 -2 -2 -1 -3 -1 -1 -2 -2 -3 -3 -1 -2 -4 7 -1 -1 -4 -3 -2 -2 -1 -2 -4+S 1 -1 1 0 -1 0 0 0 -1 -2 -2 0 -1 -2 -1 4 1 -3 -2 -2 0 0 0 -4+T 0 -1 0 -1 -1 -1 -1 -2 -2 -1 -1 -1 -1 -2 -1 1 5 -2 -2 0 -1 -1 0 -4+W -3 -3 -4 -4 -2 -2 -3 -2 -2 -3 -2 -3 -1 1 -4 -3 -2 11 2 -3 -4 -3 -2 -4+Y -2 -2 -2 -3 -2 -1 -2 -3 2 -1 -1 -2 -1 3 -3 -2 -2 2 7 -1 -3 -2 -1 -4+V 0 -3 -3 -3 -1 -2 -2 -3 -3 3 1 -2 1 -1 -2 -2 0 -3 -1 4 -3 -2 -1 -4+B -2 -1 3 4 -3 0 1 -1 0 -3 -4 0 -3 -3 -2 0 -1 -4 -3 -3 4 1 -1 -4+Z -1 0 0 1 -3 3 4 -2 0 -3 -3 1 -1 -3 -1 0 -1 -3 -2 -2 1 4 -1 -4+X 0 -1 -1 -1 -2 -1 -1 -1 -1 -1 -1 -1 -1 -1 -2 0 0 -2 -1 -1 -1 -1 -1 -4+* -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 1+|]++[matrix|PAM250+#+# This matrix was produced by "pam" Version 1.0.6 [28-Jul-93]+#+# PAM 250 substitution matrix, scale = ln(2)/3 = 0.231049+#+# Expected score = -0.844, Entropy = 0.354 bits+#+# Lowest score = -8, Highest score = 17+#+ A R N D C Q E G H I L K M F P S T W Y V B Z X *+A 2 -2 0 0 -2 0 0 1 -1 -1 -2 -1 -1 -3 1 1 1 -6 -3 0 0 0 0 -8+R -2 6 0 -1 -4 1 -1 -3 2 -2 -3 3 0 -4 0 0 -1 2 -4 -2 -1 0 -1 -8+N 0 0 2 2 -4 1 1 0 2 -2 -3 1 -2 -3 0 1 0 -4 -2 -2 2 1 0 -8+D 0 -1 2 4 -5 2 3 1 1 -2 -4 0 -3 -6 -1 0 0 -7 -4 -2 3 3 -1 -8+C -2 -4 -4 -5 12 -5 -5 -3 -3 -2 -6 -5 -5 -4 -3 0 -2 -8 0 -2 -4 -5 -3 -8+Q 0 1 1 2 -5 4 2 -1 3 -2 -2 1 -1 -5 0 -1 -1 -5 -4 -2 1 3 -1 -8+E 0 -1 1 3 -5 2 4 0 1 -2 -3 0 -2 -5 -1 0 0 -7 -4 -2 3 3 -1 -8+G 1 -3 0 1 -3 -1 0 5 -2 -3 -4 -2 -3 -5 0 1 0 -7 -5 -1 0 0 -1 -8+H -1 2 2 1 -3 3 1 -2 6 -2 -2 0 -2 -2 0 -1 -1 -3 0 -2 1 2 -1 -8+I -1 -2 -2 -2 -2 -2 -2 -3 -2 5 2 -2 2 1 -2 -1 0 -5 -1 4 -2 -2 -1 -8+L -2 -3 -3 -4 -6 -2 -3 -4 -2 2 6 -3 4 2 -3 -3 -2 -2 -1 2 -3 -3 -1 -8+K -1 3 1 0 -5 1 0 -2 0 -2 -3 5 0 -5 -1 0 0 -3 -4 -2 1 0 -1 -8+M -1 0 -2 -3 -5 -1 -2 -3 -2 2 4 0 6 0 -2 -2 -1 -4 -2 2 -2 -2 -1 -8+F -3 -4 -3 -6 -4 -5 -5 -5 -2 1 2 -5 0 9 -5 -3 -3 0 7 -1 -4 -5 -2 -8+P 1 0 0 -1 -3 0 -1 0 0 -2 -3 -1 -2 -5 6 1 0 -6 -5 -1 -1 0 -1 -8+S 1 0 1 0 0 -1 0 1 -1 -1 -3 0 -2 -3 1 2 1 -2 -3 -1 0 0 0 -8+T 1 -1 0 0 -2 -1 0 0 -1 0 -2 0 -1 -3 0 1 3 -5 -3 0 0 -1 0 -8+W -6 2 -4 -7 -8 -5 -7 -7 -3 -5 -2 -3 -4 0 -6 -2 -5 17 0 -6 -5 -6 -4 -8+Y -3 -4 -2 -4 0 -4 -4 -5 0 -1 -1 -4 -2 7 -5 -3 -3 0 10 -2 -3 -4 -2 -8+V 0 -2 -2 -2 -2 -2 -2 -1 -2 4 2 -2 2 -1 -1 -1 0 -6 -2 4 -2 -2 -1 -8+B 0 -1 2 3 -4 1 3 0 1 -2 -3 1 -2 -4 -1 0 0 -5 -3 -2 3 2 -1 -8+Z 0 0 1 3 -5 3 3 0 2 -2 -3 0 -2 -5 0 0 -1 -6 -4 -2 2 3 -1 -8+X 0 -1 0 -1 -3 -1 -1 -1 -1 -1 -1 -1 -1 -2 -1 0 0 -4 -2 -1 -1 -1 -1 -8+* -8 -8 -8 -8 -8 -8 -8 -8 -8 -8 -8 -8 -8 -8 -8 -8 -8 -8 -8 -8 -8 -8 -8 1+|]++[matrix|NUC44+#+# This matrix was created by Todd Lowe 12/10/92+#+# Uses ambiguous nucleotide codes, probabilities rounded to+# nearest integer+#+# Lowest score = -4, Highest score = 5+#+ A T G C S W R Y K M B V H D N+A 5 -4 -4 -4 -4 1 1 -4 -4 1 -4 -1 -1 -1 -2+T -4 5 -4 -4 -4 1 -4 1 1 -4 -1 -4 -1 -1 -2+G -4 -4 5 -4 1 -4 1 -4 1 -4 -1 -1 -4 -1 -2+C -4 -4 -4 5 1 -4 -4 1 -4 1 -1 -1 -1 -4 -2+S -4 -4 1 1 -1 -4 -2 -2 -2 -2 -1 -1 -3 -3 -1+W 1 1 -4 -4 -4 -1 -2 -2 -2 -2 -3 -3 -1 -1 -1+R 1 -4 1 -4 -2 -2 -1 -4 -2 -2 -3 -1 -3 -1 -1+Y -4 1 -4 1 -2 -2 -4 -1 -2 -2 -1 -3 -1 -3 -1+K -4 1 1 -4 -2 -2 -2 -2 -1 -4 -1 -3 -3 -1 -1+M 1 -4 -4 1 -2 -2 -2 -2 -4 -1 -3 -1 -1 -3 -1+B -4 -1 -1 -1 -1 -3 -3 -1 -1 -3 -1 -2 -2 -2 -1+V -1 -4 -1 -1 -1 -3 -1 -3 -3 -1 -2 -1 -2 -2 -1+H -1 -1 -4 -1 -3 -1 -3 -1 -3 -1 -2 -2 -1 -2 -1+D -1 -1 -1 -4 -3 -1 -1 -3 -1 -3 -2 -2 -2 -1 -1+N -2 -2 -2 -2 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1+|]
+ src/Bio/Chain/Alignment/Scoring/Loader.hs view
@@ -0,0 +1,30 @@+{-# LANGUAGE BangPatterns #-}+module Bio.Chain.Alignment.Scoring.Loader where++import Control.Applicative ( liftA2 )++class ScoringMatrix a where+ scoring :: a -> Char -> Char -> Int++loadMatrix :: String -> [((Char, Char), Int)]+loadMatrix txt = concatMap (lineMap . words) (tail txtlns)+ where+ -- Lines of matrix+ txtlns :: [String]+ !txtlns = filter (liftA2 (||) null ((/= '#') . head)) (strip <$> lines txt)++ -- Letters of matrix+ letters :: [Char]+ !letters = (map head . words . head) txtlns+ + -- Strip spaces+ strip :: String -> String+ strip = reverse . dropWhile (== ' ') . reverse . dropWhile (== ' ')+ + -- Map one line to a matrix+ lineMap :: [String] -> [((Char, Char), Int)]+ lineMap [] = []+ lineMap (x : xs) =+ let !hx = head x+ g (n, c) = ((hx, c), n)+ in g <$> (read <$> xs) `zip` letters
+ src/Bio/Chain/Alignment/Scoring/TH.hs view
@@ -0,0 +1,56 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE TemplateHaskell #-}+module Bio.Chain.Alignment.Scoring.TH where++import Data.Char (toLower)+import Language.Haskell.TH+import Language.Haskell.TH.Quote++import Bio.Chain.Alignment.Scoring.Loader++type Substitution a = a -> a -> Int++matrix :: QuasiQuoter+matrix = QuasiQuoter { quotePat = undefined+ , quoteType = undefined+ , quoteExp = undefined+ , quoteDec = matrixDec+ }+++matrixDec :: String -> Q [Dec]+matrixDec s = do let slines = lines s+ let txt = (unlines . tail) slines+ let name = head slines+ (typeName, dataDecl) <- typeDec name+ (funcName, funDecls) <- functionDec name txt+ smDecl <- instDec typeName funcName+ return $ dataDecl : smDecl : funDecls++instDec :: Name -> Name -> Q Dec+instDec typeN funN = return $ decl [body]+ where decl = InstanceD Nothing [] (AppT (ConT ''ScoringMatrix) (ConT typeN))+ body = FunD 'scoring [Clause [WildP] (NormalB (VarE funN)) []]++typeDec :: String -> Q (Name, Dec)+typeDec name = do typeN <- newName name+ let dataN = mkName (nameBase typeN)+#if !MIN_VERSION_template_haskell(2,12,0)+ let dervs = [ConT ''Show, ConT ''Eq]+#else+ let dervs = [DerivClause Nothing [ConT ''Show, ConT ''Eq]]+#endif+ return (typeN, DataD [] typeN [] Nothing [NormalC dataN []] dervs)++functionDec :: String -> String -> Q (Name, [Dec])+functionDec name txt = do let subM = loadMatrix txt+ funName <- newName (toLower <$> name)+ let funSign = SigD funName (AppT (ConT ''Substitution) (ConT ''Char))+ let clauses = mkClause <$> subM+ let funDecl = FunD funName clauses+ return (funName, [funSign, funDecl])++mkClause :: ((Char, Char), Int) -> Clause+mkClause ((c, d), i) = Clause [litC c, litC d] (NormalB (litI i)) []+ where litC = LitP . CharL+ litI = LitE . IntegerL . fromIntegral
+ src/Bio/Chain/Alignment/Type.hs view
@@ -0,0 +1,164 @@+{-# LANGUAGE DeriveAnyClass #-}+{-# LANGUAGE StandaloneDeriving #-}++module Bio.Chain.Alignment.Type where++import Bio.Chain (ChainLike (..))+import Control.DeepSeq (NFData (..))+import Control.Lens (Index, IxValue)+import Data.Array.Unboxed (Ix, UArray)+import GHC.Generics (Generic (..))++-- | Scoring function, returns substitution score for a couple of elements+--+-- type Scoring = Char -> Char -> Int++type Scoring a b = a -> b -> Int++-- | Simple gap penalty+--+type SimpleGap = Int++-- | Gap penalty with different 'SimpleGap' penalties for sequences.+--+-- First element of pair is penalty for first sequence passed to alignment+-- algorithm, second element — penalty for second passed sequence.+--+type SimpleGap2 = (SimpleGap, SimpleGap)++-- | Affine gap penalty+--+data AffineGap = AffineGap { gapOpen :: Int+ , gapExtend :: Int+ }+ deriving (Show, Eq, Generic, NFData)++-- | Gap penalty with different 'AffineGap' penalties for sequences.+--+-- First element of pair is penalty for first sequence passed to alignment+-- algorithm, second element — penalty for second passed sequence.+--+type AffineGap2 = (AffineGap, AffineGap)++-- | Type class that describes possible gaps in alignments.+--+class IsGap a where+ -- | Insertions are gaps in the first argument of an alignment function.+ --+ insertCostOpen :: a -> Int+ insertCostExtend :: a -> Int++ -- | Deletions are gaps in the second argument of an alignment function.+ --+ deleteCostOpen :: a -> Int+ deleteCostExtend :: a -> Int++ isAffine :: a -> Bool+ isAffine x = insertCostOpen x /= insertCostExtend x || deleteCostOpen x /= deleteCostExtend x++instance IsGap SimpleGap where+ insertCostOpen = id+ insertCostExtend = id++ deleteCostOpen = id+ deleteCostExtend = id++instance IsGap SimpleGap2 where+ insertCostOpen = fst+ insertCostExtend = fst++ deleteCostOpen = snd+ deleteCostExtend = snd++instance IsGap AffineGap where+ insertCostOpen = gapOpen+ insertCostExtend = gapExtend++ deleteCostOpen = gapOpen+ deleteCostExtend = gapExtend++instance IsGap AffineGap2 where+ insertCostOpen = gapOpen . fst+ insertCostExtend = gapExtend . fst++ deleteCostOpen = gapOpen . snd+ deleteCostExtend = gapExtend . snd++-- | Edit operation could be insertion, deletion or match/mismatch+--+data EditOp = Insert | Delete | Match+ deriving (Show, Eq, Ord, Bounded, Enum, Ix, Generic, NFData)++-- | Operation that was performed on current step of alignment+--+data Operation i j = INSERT { getJ :: j }+ | DELETE { getI :: i }+ | MATCH { getI :: i, getJ :: j }+ deriving (Show, Eq, Ord, Generic, NFData)++isInsert, isDelete, isMatch :: Operation i j -> Bool++isInsert INSERT{} = True+isInsert _ = False++isDelete DELETE{} = True+isDelete _ = False++isMatch MATCH{} = True+isMatch _ = False++-- | Alignment matrix type+--+type Matrix m m' = UArray (Index m, Index m', EditOp) Int++-- | Traceback condition type+--+type Condition m m' = Matrix m m' -> m -> m' -> Index m -> Index m' -> Bool++-- | 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?+ }++-- | Sequence Alignment result+--+data AlignmentResult m m' = AlignmentResult { score :: Int -- ^ Resulting score of alignment+ , alignment :: [Operation (Index m) (Index m')] -- ^ Alignment structure+ , sequence1 :: m -- ^ First chain+ , sequence2 :: m' -- ^ Second chain+ }+ deriving (Generic)++instance (NFData a, NFData b) => NFData (UArray (a, b, EditOp) Int) where+ rnf a = seq a ()++deriving instance (NFData a, NFData b, NFData (Index a), NFData (Index b))+ => NFData (AlignmentResult a b)++-- | Chain, that can be used for alignment+--+type Alignable m = (ChainLike m, Ix (Index m))++-- |Method of sequence alignment+--+class SequenceAlignment (a :: * -> * -> *) where+ -- | Defines wheater the alignment is semiglobal or not+ --+ semi :: a e1 e2 -> Bool+ {-# INLINABLE semi #-}+ semi = const False++ -- | Traceback conditions of alignment+ --+ cond :: (Alignable m, Alignable m') => a (IxValue m) (IxValue m') -> Conditions m m'++ -- | Starting position in matrix for traceback procedure+ --+ traceStart :: (Alignable m, Alignable m') => a (IxValue m) (IxValue m') -> Matrix m m' -> m -> m' -> (Index m, Index m')++ -- | Distance matrix element+ --+ scoreMatrix :: (Alignable m, Alignable m') => a (IxValue m) (IxValue m') -> m -> m' -> Matrix m m'
+ src/Bio/Molecule.hs view
@@ -0,0 +1,70 @@+module Bio.Molecule+ ( Molecule(..)+ , MoleculeLike(..)+ , singleton+ ) where++import Control.Lens ( (^?)+ , Index+ , IxValue+ , Ixed (..)+ , lens+ , (&)+ , (.~)+ )++newtype Molecule t c = Molecule { getChains :: [(t, c)] }+ deriving (Show, Eq)++type instance Index (Molecule t c) = t+type instance IxValue (Molecule t c) = c++class (Eq (Index m), Ixed m) => MoleculeLike m where+ -- | Create empty molecule without chains+ --+ empty :: m+ -- | Delete chain with specified index (returns error if chain doesn't present)+ --+ deleteAt :: m -> Index m -> m+ -- | Create chain with specified index (returns error if chain is already present)+ --+ create :: m -> Index m -> IxValue m -> m+ -- | Set new chain with speficied index (creates new if does not present)+ --+ set :: m -> Index m -> IxValue m -> m++-- | Create molecule with single chain+--+singleton :: MoleculeLike m => Index m -> IxValue m -> m+singleton = create empty++instance Eq t => Ixed (Molecule t c) where+ ix idx = lens (lookup idx . getChains) (\(Molecule m) my -> Molecule $ setL my m) . traverse+ where+ setL :: Maybe c -> [(t, c)] -> [(t, c)]+ setL Nothing xs = xs+ setL (Just _) [] = error "Chain should be present"+ setL y@(Just a) ((x', y') : xs) | x' == idx = (idx, a) : xs+ | otherwise = (x', y') : setL y xs++instance Eq t => MoleculeLike (Molecule t c) where+ empty = Molecule []++ deleteAt (Molecule xs) idx = Molecule $ deleteFromList xs+ where+ deleteFromList :: [(t, c)] -> [(t, c)]+ deleteFromList [] = error "Chain is not present"+ deleteFromList (a@(x', _) : ys) | x' == idx = ys+ | otherwise = a : deleteFromList ys++ create (Molecule xs) idx c = Molecule $ createInList xs+ where+ createInList :: [(t, c)] -> [(t, c)]+ createInList [] = [(idx, c)]+ createInList (a@(x', _) : ys)+ | x' == idx = error "Chain should not be present at molecule"+ | otherwise = a : createInList ys++ set m idx c = case m ^? ix idx of+ Nothing -> create m idx c+ Just _ -> m & ix idx .~ c
+ src/Bio/NucleicAcid/Chain.hs view
@@ -0,0 +1,24 @@+{-# OPTIONS_GHC -fno-warn-orphans #-}++module Bio.NucleicAcid.Chain where++import Bio.Chain as C+import Bio.NucleicAcid.Nucleotide+import Control.Lens+import Data.Array (Ix (..))+import Data.String (IsString (..))++newtype NucleicAcidChain i a = NucleicAcidChain { getChain :: Chain i a }+ deriving (Show, Eq, Functor, Foldable, Traversable, ChainLike)++type instance Index (NucleicAcidChain i a) = i+type instance IxValue (NucleicAcidChain i a) = a++instance Ix i => Ixed (NucleicAcidChain i a) where+ ix i' = coerced . ix @(Chain i a) i'++instance IsString (NucleicAcidChain Int DNA) where+ fromString = NucleicAcidChain . fromString++instance IsString (NucleicAcidChain Int RNA) where+ fromString = NucleicAcidChain . fromString
+ src/Bio/NucleicAcid/Nucleotide.hs view
@@ -0,0 +1,8 @@+module Bio.NucleicAcid.Nucleotide+ ( module T+ , module M+ ) where++import Bio.NucleicAcid.Nucleotide.Instances ()+import Bio.NucleicAcid.Nucleotide.Type as T+import Bio.Utils.Monomer as M
+ src/Bio/NucleicAcid/Nucleotide/Instances.hs view
@@ -0,0 +1,60 @@+{-# OPTIONS_GHC -fno-warn-orphans #-}++module Bio.NucleicAcid.Nucleotide.Instances () where++import Bio.NucleicAcid.Nucleotide.Type+import Bio.Utils.Monomer (FromSymbol (..), Symbol (..))+import Data.Array (Array, listArray)+import Data.String (IsString (..))++-------------------------------------------------------------------------------+-- Symbol and ThreeSymbols+-------------------------------------------------------------------------------++instance Symbol DNA where+ symbol DA = 'A'+ symbol DC = 'C'+ symbol DG = 'G'+ symbol DT = 'T'++instance FromSymbol DNA where+ fromSymbolE 'A' = Right DA+ fromSymbolE 'C' = Right DC+ fromSymbolE 'G' = Right DG+ fromSymbolE 'T' = Right DT+ fromSymbolE ch = Left ch++instance Symbol RNA where+ symbol RA = 'A'+ symbol RC = 'C'+ symbol RG = 'G'+ symbol RU = 'U'++instance FromSymbol RNA where+ fromSymbolE 'A' = Right RA+ fromSymbolE 'C' = Right RC+ fromSymbolE 'G' = Right RG+ fromSymbolE 'U' = Right RU+ fromSymbolE ch = Left ch++-------------------------------------------------------------------------------+-- IsString+-------------------------------------------------------------------------------++instance {-# OVERLAPPING #-} IsString [DNA] where+ fromString s =+ case traverse fromSymbolE s of+ Right l -> l+ Left e -> error $ "Bio.NucleicAcid.Nucleotide.Instances: could not read nucleotide " <> [e]++instance IsString (Array Int DNA) where+ fromString s = listArray (0, length s - 1) $ fromString s++instance {-# OVERLAPPING #-} IsString [RNA] where+ fromString s =+ case traverse fromSymbolE s of+ Right l -> l+ Left e -> error $ "Bio.NucleicAcid.Nucleotide.Instances: could not read nucleotide " <> [e]++instance IsString (Array Int RNA) where+ fromString s = listArray (0, length s - 1) $ fromString s
+ src/Bio/NucleicAcid/Nucleotide/Type.hs view
@@ -0,0 +1,94 @@+{-# LANGUAGE DeriveAnyClass #-}++module Bio.NucleicAcid.Nucleotide.Type+ ( DNA (..)+ , RNA (..)+ , nucleoIso+ , toRNA+ , toDNA+ , Complementary (..)+ ) where++import Control.DeepSeq (NFData)+import Control.Lens (Iso', iso)+import Data.Array (Array, Ix, bounds, listArray)+import Data.Foldable (Foldable (..))+import GHC.Generics (Generic)++data DNA = DA | DC | DG | DT+ deriving (Eq, Ord, Bounded, Enum, Generic, NFData)++instance Show DNA where+ show DA = "Adenine"+ show DC = "Cytosine"+ show DG = "Guanine"+ show DT = "Thymine"++data RNA = RA | RC | RG | RU+ deriving (Eq, Ord, Bounded, Enum, Generic, NFData)++instance Show RNA where+ show RA = "Adenine"+ show RC = "Cytosine"+ show RG = "Guanine"+ show RU = "Uracil"++-------------------------------------------------------------------------------+-- Transciption+-------------------------------------------------------------------------------++nucleoIso :: Iso' DNA RNA+nucleoIso = iso toRNA toDNA++{-# INLINE toRNA #-}+toRNA :: DNA -> RNA+toRNA DA = RA+toRNA DC = RC+toRNA DG = RG+toRNA DT = RU++{-# INLINE toDNA #-}+toDNA :: RNA -> DNA+toDNA RA = DA+toDNA RC = DC+toDNA RG = DG+toDNA RU = DT++------------------------------------------------------------------------------+-- Complementary and reverse complementary+-------------------------------------------------------------------------------++class Complementary a where+ -- | complement *NA (DNA or RNA)+ --+ cNA :: a -> a++ -- | reverce complement *NA (DNA or RNA)+ --+ rcNA :: a -> a++instance Complementary DNA where+ cNA DA = DT+ cNA DC = DG+ cNA DG = DC+ cNA DT = DA++ rcNA = cNA++instance Complementary RNA where+ cNA = toRNA . cNA . toDNA++ rcNA = cNA++instance Complementary a => Complementary [a] where+ cNA = fmap cNA++ rcNA = reverse . cNA++instance (Complementary a, Ix i) => Complementary (Array i a) where+ cNA = fmap cNA++ rcNA l = listArray (bounds l) rl+ where+ rl = rcNA . toList $ l+
+ src/Bio/Protein/Algebra.hs view
@@ -0,0 +1,203 @@+module Bio.Protein.Algebra+ ( phi+ , psi+ , omega+ , chi+ ) where++import Data.Monoid ( First (..) )+import Control.Lens+import Bio.Utils.Geometry ( V3R, R, Ray (..), normalize, rotateR )++import Bio.Protein.AminoAcid+import Bio.Protein.Metric+import Bio.Protein.Chain++type Dihedral m f r g h = (ChainLike m, HasN f, HasCA r, HasC g, HasAtom h, IxValue m ~ AminoAcid f r g (h V3R))++-- | Measure and rotate Psi dihedral angle+--+psi :: forall m f r g h.Dihedral m f r g h => Index m -> Traversal' m R+psi i = rcd (\rot -> (& c %~ fmap rot)) (ix i . n . atom) (ix i . ca . atom) (ix i . c . atom) (ix (succ i) . n . atom) i++-- | Measure and rotate Phi dihedral angle+--+phi :: forall m f r g h.Dihedral m f r g h => Index m -> Traversal' m R+phi i = rcd (\rot -> (& ca %~ fmap rot) . (& c %~ fmap rot)) (ix (pred i) . c . atom) (ix i . n . atom) (ix i . ca . atom) (ix i . c . atom) i++-- | Measure and rotate Omega dihedral angle+--+omega :: forall m f r g h.Dihedral m f r g h => Index m -> Traversal' m R+omega i = rcd (fmap . fmap) (ix (pred i) . ca . atom) (ix (pred i) . c . atom) (ix i . n . atom) (ix i . ca . atom) i++-- | Measure and rotate Chi (1, 2, 3, 4, 5) dihedral angles+--+chi :: forall nr cr h m.(HasN nr, Functor cr, HasAtom h, m ~ AminoAcid nr (Env Radical) cr (h V3R)) => Int -> Traversal' m R+chi i = lens getChi setChi . traverse+ where+ checkI :: Bool+ checkI = i > 0 && i < 6++ getChi :: m -> Maybe R+ getChi | checkI = (^? dihedral @(First V3R) (chiP i) (chiP (i + 1)) (chiP (i + 2)) (chiP (i + 3)))+ | otherwise = const Nothing++ setChi :: m -> Maybe R -> m+ setChi m Nothing = m+ setChi m (Just d) | checkI = safeSetChi m d+ | otherwise = m++ safeSetChi :: m -> R -> m+ safeSetChi m d = case getChi m of+ Nothing -> m+ Just cud ->+ let ray = Ray (m ^?! chiP (i + 1)) (normalize $ m ^?! chiP (i + 2) - m ^?! chiP (i + 1))+ rot = rotateR ray (cud - d) :: V3R -> V3R+ in rotateRadical i rot m++ rotateRadical :: Int -> (V3R -> V3R) -> m -> m+ rotateRadical j rot m | j == 1 && m ^. radicalType /= PRO = rr $ m & radical %~ fmap (fmap rot)+ -- Chi 2+ | j == 2 && m ^. radicalType == ASP = m & radical . cg %~ fmap rot+ & radical . od1 %~ fmap rot+ & radical . od2 %~ fmap rot+ | j == 2 && m ^. radicalType == PHE = m & radical . cg %~ fmap rot+ & radical . cd1 %~ fmap rot+ & radical . cd2 %~ fmap rot+ & radical . ce1 %~ fmap rot+ & radical . ce2 %~ fmap rot+ & radical . cz %~ fmap rot+ | j == 2 && m ^. radicalType == HIS = m & radical . cg %~ fmap rot+ & radical . nd1 %~ fmap rot+ & radical . cd2 %~ fmap rot+ & radical . ce1 %~ fmap rot+ & radical . ne2 %~ fmap rot+ | j == 2 && m ^. radicalType == ILE = m & radical . cg1 %~ fmap rot+ & radical . cg2 %~ fmap rot+ & radical . cd1 %~ fmap rot+ | j == 2 && m ^. radicalType == LEU = m & radical . cg %~ fmap rot+ & radical . cd1 %~ fmap rot+ & radical . cd2 %~ fmap rot+ | j == 2 && m ^. radicalType == ASN = m & radical . cg %~ fmap rot+ & radical . od1 %~ fmap rot+ & radical . nd2 %~ fmap rot+ | j == 2 && m ^. radicalType == TRP = m & radical . cg %~ fmap rot+ & radical . cd1 %~ fmap rot+ & radical . cd2 %~ fmap rot+ & radical . ne1 %~ fmap rot+ & radical . ce2 %~ fmap rot+ & radical . ce3 %~ fmap rot+ & radical . cz2 %~ fmap rot+ & radical . cz3 %~ fmap rot+ & radical . ch2 %~ fmap rot+ | j == 2 && m ^. radicalType == TYR = m & radical . cg %~ fmap rot+ & radical . cd1 %~ fmap rot+ & radical . cd2 %~ fmap rot+ & radical . ce1 %~ fmap rot+ & radical . ce2 %~ fmap rot+ & radical . cz %~ fmap rot+ & radical . ch2 %~ fmap rot+ | j == 2 && m ^. radicalType == GLU = rr $ m & radical . cg %~ fmap rot+ | j == 2 && m ^. radicalType == MET = rr $ m & radical . cg %~ fmap rot+ | j == 2 && m ^. radicalType == GLN = rr $ m & radical . cg %~ fmap rot+ | j == 2 && m ^. radicalType == LYS = rr $ m & radical . cg %~ fmap rot+ | j == 2 && m ^. radicalType == ARG = rr $ m & radical . cg %~ fmap rot+ -- Chi 3+ | j == 3 && m ^. radicalType == GLU = m & radical . cd %~ fmap rot+ & radical . oe1 %~ fmap rot+ & radical . oe2 %~ fmap rot+ | j == 3 && m ^. radicalType == MET = m & radical . sd %~ fmap rot+ & radical . ce %~ fmap rot+ | j == 3 && m ^. radicalType == GLN = m & radical . cd %~ fmap rot+ & radical . oe1 %~ fmap rot+ & radical . ne2 %~ fmap rot+ | j == 3 && m ^. radicalType == LYS = rr $ m & radical . cd %~ fmap rot+ | j == 3 && m ^. radicalType == ARG = rr $ m & radical . cd %~ fmap rot+ -- Chi 4+ | j == 4 && m ^. radicalType == LYS = m & radical . ce %~ fmap rot+ & radical . nz %~ fmap rot+ | j == 4 && m ^. radicalType == ARG = rr $ m & radical . ne %~ fmap rot+ -- Chi 5+ | j == 5 && m ^. radicalType == ARG = m & radical . cz %~ fmap rot+ & radical . nh1 %~ fmap rot+ & radical . nh2 %~ fmap rot+ | otherwise = m+ where+ rr = rotateRadical (j + 1) rot++-- Helper functions++-- | Chi angle point (one to eight)+-- Points 1, 2, 3 and 4 — Chi 1+-- Points 2, 3, 4 and 5 - Chi 2+-- Points 3, 4, 5 and 6 - Chi 3+-- Points 4, 5, 6 and 7 - Chi 4+-- Points 5, 6, 7 and 8 - Chi 5+--+chiP :: forall a nr cr h m.(HasN nr, Functor cr, HasAtom h, m ~ AminoAcid nr (Env Radical) cr (h a)) => Int -> Traversal' m a+chiP i = lens getChiP setChiP . traverse+ where+ checkI :: Bool+ checkI = i > 0 && i < 9++ chiPL :: Int -> AA -> Traversal' m a+ chiPL 1 _ = n . atom+ chiPL 2 _ = ca . atom+ chiPL 3 _ = radical . cb . atom+ chiPL 4 aa | aa == CYS = radical . sg . atom+ | aa == ILE = radical . cg1 . atom+ | aa == SER = radical . og . atom+ | aa == THR = radical . og1 . atom+ | otherwise = radical . cg . atom+ chiPL 5 aa | aa == ASN = radical . od1 . atom+ | aa == ASP = radical . od1 . atom+ | aa == HIS = radical . nd1 . atom+ | aa == MET = radical . sd . atom+ | aa == LEU = radical . cd1 . atom+ | aa == PHE = radical . cd1 . atom+ | aa == TRP = radical . cd1 . atom+ | aa == TYR = radical . cd1 . atom+ | otherwise = radical . cd . atom+ chiPL 6 aa | aa == ARG = radical . ne . atom+ | aa == GLN = radical . oe1 . atom+ | aa == GLU = radical . oe1 . atom+ | otherwise = radical . ce . atom+ chiPL 7 aa | aa == LYS = radical . nz . atom+ | otherwise = radical . cz . atom+ chiPL 8 _ = radical . nh1 . atom+ chiPL _ _ = error "You cannot be here, as Chi dihedrals involves only 8 points"++ getChiP :: m -> Maybe a+ getChiP m | checkI = m ^? chiPL i (m ^. radicalType)+ | otherwise = Nothing++ setChiP :: m -> Maybe a -> m+ setChiP m Nothing = m+ setChiP m (Just v) | checkI = over (chiPL i (m ^. radicalType)) (const v) m+ | otherwise = m++type ModifyFunction m = (V3R -> V3R) -> IxValue m -> IxValue m++-- | Rotate cannonical dihedral in backbone+--+rcd :: forall m f r g h.Dihedral m f r g h => ModifyFunction m {- modify function -} ->+ Traversal' m V3R {- first point -} ->+ Traversal' m V3R {- second point -} ->+ Traversal' m V3R {- third point -} ->+ Traversal' m V3R {- fourth point -} ->+ Index m {- dihedral index -} ->+ Traversal' m R+rcd mf x1 x2 x3 x4 i = lens getRCD setRCD . traverse+ where+ getRCD :: m -> Maybe R+ getRCD = (^? dihedral @(First V3R) x1 x2 x3 x4)++ setRCD :: m -> Maybe R -> m+ setRCD ar Nothing = ar+ setRCD ar (Just d) = case getRCD ar of+ Nothing -> ar+ Just cud ->+ let ray = Ray (ar ^?! x2) (normalize $ ar ^?! x3 - ar ^?! x2)+ rot = rotateR ray (cud - d)+ mfy = modify i (mf rot) . modifyAfter i (fmap (fmap rot))+ in mfy ar
+ src/Bio/Protein/AminoAcid.hs view
@@ -0,0 +1,9 @@+module Bio.Protein.AminoAcid+ ( module T+ , module I+ , module M+ ) where++import Bio.Protein.AminoAcid.Type as T+import Bio.Protein.AminoAcid.Instances as I+import Bio.Utils.Monomer as M
+ src/Bio/Protein/AminoAcid/Instances.hs view
@@ -0,0 +1,326 @@+{-# OPTIONS_GHC -fno-warn-orphans #-}+{-# LANGUAGE OverloadedStrings #-}++module Bio.Protein.AminoAcid.Instances where++import Bio.Protein.AminoAcid.Type+import Bio.Utils.Monomer (FromSymbol (..),+ FromThreeSymbols (..), Symbol (..),+ ThreeSymbols (..))+import Control.Lens (Const (..), Getting, Identity (..),+ Lens', coerced, to, (^.))+import Data.Array (Array, listArray)+import Data.Coerce (coerce)+import Data.String (IsString (..))++-------------------------------------------------------------------------------+-- Creatable+-------------------------------------------------------------------------------++-- | Single object can be created+--+class Createable a where+ type Create a :: *+ -- | Function to create single object+ --+ create :: Create a++instance Createable (BB a) where+ type Create (BB a) = a -> a -> a -> BB a+ create n_ ca_ c_ = coerce <$> AminoAcid (pure n_) (pure ca_) (pure c_)++instance Createable (BBCA a) where+ type Create (BBCA a) = a -> BBCA a+ create ca_ = coerce <$> AminoAcid (Const ()) (pure ca_) (Const ())++instance Createable (BBT a) where+ type Create (BBT a) = a -> a -> a -> AA -> BBT a+ create n_ ca_ c_ aa = coerce <$> AminoAcid (pure n_) (Env ca_ (Const aa)) (pure c_)++instance Createable (BBCAT a) where+ type Create (BBCAT a) = a -> AA -> BBCAT a+ create ca_ aa = coerce <$> AminoAcid (Const ()) (Env ca_ (Const aa)) (Const ())++instance Createable (BBCG a) where+ type Create (BBCG a) = a -> a -> a -> a -> AA -> BBCG a+ create n_ ca_ c_ cg_ aa = coerce <$> AminoAcid (pure n_) (Env ca_ (CG cg_ aa)) (pure c_)++instance Createable (BBO a) where+ type Create (BBO a) = a -> a -> a -> a -> BBO a+ create n_ ca_ c_ o_ = coerce <$> AminoAcid (pure n_) (pure ca_) (Env c_ (pure o_))++instance Createable (BBOT a) where+ type Create (BBOT a) = a -> a -> a -> a -> AA -> BBOT a+ create n_ ca_ c_ o_ aa = coerce <$> AminoAcid (pure n_) (Env ca_ (Const aa)) (Env c_ (pure o_))++instance Createable (BBOCG a) where+ type Create (BBOCG a) = a -> a -> a -> a -> a -> AA -> BBOCG a+ create n_ ca_ c_ o_ cg_ aa = coerce <$> AminoAcid (pure n_) (Env ca_ (CG cg_ aa)) (Env c_ (pure o_))++instance Createable (BBOR a) where+ type Create (BBOR a) = a -> a -> a -> a -> Radical a -> BBOR a+ create n_ ca_ c_ o_ r = coerce <$> AminoAcid (pure n_) (Env ca_ r) (Env c_ (pure o_))++instance Createable (BBOXTR a) where+ type Create (BBOXTR a) = a -> a -> a -> a -> a -> Radical a -> BBOXTR a+ create n_ ca_ c_ o_ oxt_ r = coerce <$> AminoAcid (pure n_) (Env ca_ r) (Env c_ (OXT o_ oxt_))++instance Createable (BBORH a) where+ type Create (BBORH a) = a -> a -> a -> a -> Radical a -> BBORH a+ create n_ ca_ c_ o_ r = flip Env [] <$> AminoAcid (pure n_) (Env ca_ r) (Env c_ (pure o_))++instance Createable (BBOXTRH a) where+ type Create (BBOXTRH a) = a -> a -> a -> a -> a -> Radical a -> BBOXTRH a+ create n_ ca_ c_ o_ oxt_ r = flip Env [] <$> AminoAcid (pure n_) (Env ca_ r) (Env c_ (OXT o_ oxt_))++-------------------------------------------------------------------------------+-- HasRadical+-------------------------------------------------------------------------------++-- | Has lens to observe, set and modify radicals+--+class Functor r => HasRadical r where+ type RadicalType r a :: *+ -- | Lens for radical atom or group+ --+ radical :: (Functor f, Functor g) => Lens' (AminoAcid f (Env r) g a) (RadicalType r a)++instance HasRadical (Const x) where+ type RadicalType (Const x) a = x+ radical = ca' . environment . coerced++instance HasRadical Radical where+ type RadicalType Radical a = Radical a+ radical = ca' . environment++instance HasRadical CG where+ type RadicalType CG a = a+ radical = ca' . environment . cg'++instance HasRadical Identity where+ type RadicalType Identity a = a+ radical = ca' . environment . coerced++-------------------------------------------------------------------------------+-- HasRadicalType+-------------------------------------------------------------------------------++-- | Has lens to observe radical types+--+class Functor r => HasRadicalType r where+ -- | Getter for radical type+ --+ radicalType :: (Functor f, Functor g) => Getting AA (AminoAcid f (Env r) g a) AA++instance HasRadicalType (Const AA) where+ radicalType = ca' . environment . coerced++instance HasRadicalType CG where+ radicalType = ca' . environment . radical' . coerced++instance HasRadicalType Radical where+ radicalType = ca' . environment . to rad2rad++-------------------------------------------------------------------------------+-- Has some atom+-------------------------------------------------------------------------------++-- | Has lens to observe, set and modify ca_ atom+--+class Functor r => HasCA r where+ -- | Lens for ca_ atom+ --+ ca :: (Functor f, Functor g) => Lens' (AminoAcid f r g a) a++instance HasCA Identity where+ ca = ca' . coerced++instance Functor f => HasCA (Env f) where+ ca = ca' . atom'++-- | Has lens to observe, set and modify c_ atom+--+class Functor r => HasC r where+ -- | Lens for c_ atom+ --+ c :: (Functor f, Functor g) => Lens' (AminoAcid f g r a) a++instance HasC Identity where+ c = c' . coerced++instance Functor f => HasC (Env f) where+ c = c' . atom'++-- | Has lens to observe, set and modify o_ atom+--+class Functor r => HasO r where+ -- | Lens for o_ atom+ --+ o :: (Functor f, Functor g) => Lens' (AminoAcid f g (Env r) a) a++instance HasO Identity where+ o = c' . environment . coerced++instance HasO OXT where+ o = c' . environment . o'++-- | Has lens to observe, set and modify OXT atom+--+class Functor r => HasOXT r where+ -- | Lens for OXT atom+ --+ oxt :: (Functor f, Functor g) => Lens' (AminoAcid f g (Env r) a) a++instance HasOXT OXT where+ oxt = c' . environment . oxt'++-- | Has lens to observe, set and modify n_ atom+--+class Functor r => HasN r where+ -- | Lens for n_ atom+ --+ n :: (Functor f, Functor g) => Lens' (AminoAcid r f g a) a++instance HasN Identity where+ n = n' . coerced++instance Functor f => HasN (Env f) where+ n = n' . atom'++-- | Lens to get atom from some enviroment+--+class Functor f => HasAtom f where+ -- | Lens for exact atom get+ --+ atom :: Lens' (f a) a++instance HasAtom Identity where+ atom = coerced++instance Functor r => HasAtom (Env r) where+ atom = atom'++-- | Lens to get hydrogens from hydrated atom+--+hydrogens :: Lens' (Env [] a) [a]+hydrogens = environment++-------------------------------------------------------------------------------+-- Symbol and ThreeSymbols+-------------------------------------------------------------------------------++-- | Lens to get Symbol from every suitable amino acid+--+instance (Functor nr, HasRadicalType car, Functor cr) => Symbol (AminoAcid nr (Env car) cr a) where+ symbol = symbol . (^. radicalType)++-- | Symbol encoding+--+instance Symbol AA where+ symbol ALA = 'A'+ symbol CYS = 'C'+ symbol ASP = 'D'+ symbol GLU = 'E'+ symbol PHE = 'F'+ symbol GLY = 'G'+ symbol HIS = 'H'+ symbol ILE = 'I'+ symbol LYS = 'K'+ symbol LEU = 'L'+ symbol MET = 'M'+ symbol ASN = 'N'+ symbol PRO = 'P'+ symbol GLN = 'Q'+ symbol ARG = 'R'+ symbol SER = 'S'+ symbol THR = 'T'+ symbol VAL = 'V'+ symbol TRP = 'W'+ symbol TYR = 'Y'++-- | Parse symbol encoding+--+instance FromSymbol AA where+ fromSymbolE 'A' = Right ALA+ fromSymbolE 'C' = Right CYS+ fromSymbolE 'D' = Right ASP+ fromSymbolE 'E' = Right GLU+ fromSymbolE 'F' = Right PHE+ fromSymbolE 'G' = Right GLY+ fromSymbolE 'H' = Right HIS+ fromSymbolE 'I' = Right ILE+ fromSymbolE 'K' = Right LYS+ fromSymbolE 'L' = Right LEU+ fromSymbolE 'M' = Right MET+ fromSymbolE 'N' = Right ASN+ fromSymbolE 'P' = Right PRO+ fromSymbolE 'Q' = Right GLN+ fromSymbolE 'R' = Right ARG+ fromSymbolE 'S' = Right SER+ fromSymbolE 'T' = Right THR+ fromSymbolE 'V' = Right VAL+ fromSymbolE 'W' = Right TRP+ fromSymbolE 'Y' = Right TYR+ fromSymbolE ch = Left ch++-- | Three symbols encoding+--+instance ThreeSymbols AA where+ threeSymbols ALA = "ALA"+ threeSymbols CYS = "CYS"+ threeSymbols ASP = "ASP"+ threeSymbols GLU = "GLU"+ threeSymbols PHE = "PHE"+ threeSymbols GLY = "GLY"+ threeSymbols HIS = "HIS"+ threeSymbols ILE = "ILE"+ threeSymbols LYS = "LYS"+ threeSymbols LEU = "LEU"+ threeSymbols MET = "MET"+ threeSymbols ASN = "ASN"+ threeSymbols PRO = "PRO"+ threeSymbols GLN = "GLN"+ threeSymbols ARG = "ARG"+ threeSymbols SER = "SER"+ threeSymbols THR = "THR"+ threeSymbols VAL = "VAL"+ threeSymbols TRP = "TRP"+ threeSymbols TYR = "TYR"++-- | Parse three symbols encoding+--+instance FromThreeSymbols AA where+ fromThreeSymbols "ALA" = Just ALA+ fromThreeSymbols "CYS" = Just CYS+ fromThreeSymbols "ASP" = Just ASP+ fromThreeSymbols "GLU" = Just GLU+ fromThreeSymbols "PHE" = Just PHE+ fromThreeSymbols "GLY" = Just GLY+ fromThreeSymbols "HIS" = Just HIS+ fromThreeSymbols "ILE" = Just ILE+ fromThreeSymbols "LYS" = Just LYS+ fromThreeSymbols "LEU" = Just LEU+ fromThreeSymbols "MET" = Just MET+ fromThreeSymbols "ASN" = Just ASN+ fromThreeSymbols "PRO" = Just PRO+ fromThreeSymbols "GLN" = Just GLN+ fromThreeSymbols "ARG" = Just ARG+ fromThreeSymbols "SER" = Just SER+ fromThreeSymbols "THR" = Just THR+ fromThreeSymbols "VAL" = Just VAL+ fromThreeSymbols "TRP" = Just TRP+ fromThreeSymbols "TYR" = Just TYR+ fromThreeSymbols _ = Nothing++-------------------------------------------------------------------------------+-- IsString+-------------------------------------------------------------------------------++instance {-# OVERLAPPING #-} IsString [AA] where+ fromString s =+ case traverse fromSymbolE s of+ Right l -> l+ Left e -> error $ "Bio.Protein.AminoAcid.Instances: could not read aminoacid " <> [e]++instance IsString (Array Int AA) where+ fromString s = listArray (0, length s - 1) $ fromString s
+ src/Bio/Protein/AminoAcid/Type.hs view
@@ -0,0 +1,299 @@+{-# LANGUAGE DeriveFunctor #-}+{-# LANGUAGE TemplateHaskell #-}+{-# LANGUAGE DeriveAnyClass #-}++module Bio.Protein.AminoAcid.Type where++import Control.DeepSeq (NFData (..))+import Control.Lens+import Control.Monad.Identity (Identity)+import GHC.Generics (Generic (..))++-- | Proteinogenic amino acids+--+data AA = ALA -- A+ | CYS -- C+ | ASP -- D+ | GLU -- E+ | PHE -- F+ | GLY -- G+ | HIS -- H+ | ILE -- I+ | LYS -- K+ | LEU -- L+ | MET -- M+ | ASN -- N+ | PRO -- P+ | GLN -- Q+ | ARG -- R+ | SER -- S+ | THR -- T+ | VAL -- V+ | TRP -- W+ | TYR -- Y+ deriving (Eq, Ord, Bounded, Enum, Generic, NFData)++-- | Show full names of amino acids+--+instance Show AA where+ show ALA = "Alanine"+ show CYS = "Cysteine"+ show ASP = "AsparticAcid"+ show GLU = "GlutamicAcid"+ show PHE = "Phenylalanine"+ show GLY = "Glycine"+ show HIS = "Histidine"+ show ILE = "Isoleucine"+ show LYS = "Lysine"+ show LEU = "Leucine"+ show MET = "Methionine"+ show ASN = "Asparagine"+ show PRO = "Proline"+ show GLN = "Glutamine"+ show ARG = "Arginine"+ show SER = "Serine"+ show THR = "Threonine"+ show VAL = "Valine"+ show TRP = "Tryptophan"+ show TYR = "Tyrosine"++-- | Amino acid structure type+--+data AminoAcid nr car cr a = AminoAcid { _n' :: nr a+ , _ca' :: car a+ , _c' :: cr a+ }+ deriving (Show, Eq, Functor, Generic, NFData)++-- | Radical structure type+--+data Radical a = Alanine -- no chi+ { _cb :: a -- -CB+ } --+ | Cysteine --+ { _cb :: a -- -CB-SG+ , _sg :: a --+ } --+ | AsparticAcid --+ { _cb :: a -- -CB-CG-OD1+ , _cg :: a -- |+ , _od1 :: a -- OD2+ , _od2 :: a --+ } --+ | GlutamicAcid --+ { _cb :: a -- -CB-CG-CD-OE1+ , _cg :: a -- |+ , _cd :: a -- OE2+ , _oe1 :: a --+ , _oe2 :: a --+ } --+ | Phenylalanine --+ { _cb :: a -- -CB-CG-CD1-CE1+ , _cg :: a -- | |+ , _cd1 :: a -- CD2-CE2-CZ+ , _cd2 :: a --+ , _ce1 :: a --+ , _ce2 :: a --+ , _cz :: a --+ } --+ | Glycine --+ | Histidine --+ { _cb :: a -- -CB-CG-ND1-CE1+ , _cg :: a -- |+ , _nd1 :: a -- CD2-NE2+ , _cd2 :: a --+ , _ce1 :: a --+ , _ne2 :: a --+ } --+ | Isoleucine --+ { _cb :: a -- -CB-CG1-CD1+ , _cg1 :: a -- |+ , _cg2 :: a -- CG2+ , _cd1 :: a --+ } --+ | Lysine --+ { _cb :: a -- -CB-CG-CD-CE-NZ+ , _cg :: a --+ , _cd :: a --+ , _ce :: a --+ , _nz :: a --+ } --+ | Leucine --+ { _cb :: a -- -CB-CG-CD1+ , _cg :: a -- |+ , _cd1 :: a -- CD2+ , _cd2 :: a --+ } --+ | Methionine --+ { _cb :: a -- -CB-CG-SD-CE+ , _cg :: a --+ , _sd :: a --+ , _ce :: a --+ } --+ | Asparagine --+ { _cb :: a -- -CB-CG-OD1+ , _cg :: a -- |+ , _od1 :: a -- ND2+ , _nd2 :: a --+ } --+ | Proline --+ { _cb :: a -- -CB-CG-CD(-N)+ , _cg :: a --+ , _cd :: a --+ } --+ | Glutamine --+ { _cb :: a -- -CB-CG-CD-OE1+ , _cg :: a -- |+ , _cd :: a -- NE2+ , _oe1 :: a --+ , _ne2 :: a --+ } --+ | Arginine --+ { _cb :: a -- -CB-CG-CD-NE-CZ-NH1+ , _cg :: a -- |+ , _cd :: a -- NH2+ , _ne :: a --+ , _cz :: a --+ , _nh1 :: a --+ , _nh2 :: a --+ } --+ | Serine --+ { _cb :: a -- -CB-OG+ , _og :: a --+ } --+ | Threonine --+ { _cb :: a -- -CB-OG1+ , _og1 :: a -- |+ , _cg2 :: a -- CG2+ } --+ | Valine --+ { _cb :: a -- -CB-CG1+ , _cg1 :: a -- |+ , _cg2 :: a -- CG2+ } --+ | Tryptophan --+ { _cb :: a -- -CB-CG-CD1-NE1+ , _cg :: a -- | |+ , _cd1 :: a -- CD2-CE2-CZ2+ , _cd2 :: a -- | |+ , _ne1 :: a -- CE3-CZ3-CH2+ , _ce2 :: a --+ , _ce3 :: a --+ , _cz2 :: a --+ , _cz3 :: a --+ , _ch2 :: a --+ } --+ | Tyrosine --+ { _cb :: a -- -CB-CG-CD1-CE1+ , _cg :: a -- | |+ , _cd1 :: a -- CD2-CE2-CZ-OH+ , _cd2 :: a --+ , _ce1 :: a --+ , _ce2 :: a --+ , _cz :: a --+ , _oh :: a --+ } --+ deriving (Show, Eq, Functor, Generic, NFData)++-- | Atom environment, e.g. hydrogens or radicals+--+data Env r a = Env { _atom' :: a+ , _environment :: r a+ }+ deriving (Show, Eq, Functor, Generic, NFData)++-- | Hydrogens envrironment+--+type H a = Env [] a++-- | Oxigen and hydroxi group, connected to C-terminal of amino acid+--+data OXT a = OXT { _o' :: a+ , _oxt' :: a+ }+ deriving (Show, Eq, Functor, Generic, NFData)++-- | CG atom with radical type+--+data CG a = CG { _cg' :: a+ , _radical' :: AA+ }+ deriving (Show, Eq, Functor, Generic, NFData)++makeLenses ''AminoAcid+makeLenses ''Radical+makeLenses ''Env+makeLenses ''OXT+makeLenses ''CG++-- | BackBone+--+type BB a = AminoAcid Identity Identity Identity (Identity a)++-- | BackBone CA-only+--+type BBCA a = AminoAcid (Const ()) Identity (Const ()) (Identity a)++-- | BackBone with radical Type+--+type BBT a = AminoAcid Identity (Env (Const AA)) Identity (Identity a)++-- | BackBone CA-only with radical Type+--+type BBCAT a = AminoAcid (Const ()) (Env (Const AA)) (Const ()) (Identity a)++-- | BackBone with CG-radical+--+type BBCG a = AminoAcid Identity (Env CG) Identity (Identity a)++-- | BackBone with Oxigen+--+type BBO a = AminoAcid Identity Identity (Env Identity) (Identity a)++-- | BackBone with Oxigen and radical Type+--+type BBOT a = AminoAcid Identity (Env (Const AA)) (Env Identity) (Identity a)++-- | BackBone with Oxigen and CG-radical+--+type BBOCG a = AminoAcid Identity (Env CG) (Env Identity) (Identity a)++-- | BackBone with Oxigen and Radical+--+type BBOR a = AminoAcid Identity (Env Radical) (Env Identity) (Identity a)++-- | BackBone with Oxigen, oXigen Two and Radical+--+type BBOXTR a = AminoAcid Identity (Env Radical) (Env OXT) (Identity a)++-- | BackBone with Oxigen, Radical and Hydrogens+--+type BBORH a = AminoAcid Identity (Env Radical) (Env Identity) (H a)++-- | BackBone with Oxigen, oXigen Two, Radical and Hydrogens+--+type BBOXTRH a = AminoAcid Identity (Env Radical) (Env OXT) (H a)+--+-- | Convert radical to radical name+--+rad2rad :: Radical a -> AA+rad2rad Alanine{} = ALA+rad2rad Cysteine{} = CYS+rad2rad AsparticAcid{} = ASP+rad2rad GlutamicAcid{} = GLU+rad2rad Phenylalanine{} = PHE+rad2rad Glycine = GLY+rad2rad Histidine{} = HIS+rad2rad Isoleucine{} = ILE+rad2rad Lysine{} = LYS+rad2rad Leucine{} = LEU+rad2rad Methionine{} = MET+rad2rad Asparagine{} = ASN+rad2rad Proline{} = PRO+rad2rad Glutamine{} = GLN+rad2rad Arginine{} = ARG+rad2rad Serine{} = SER+rad2rad Threonine{} = THR+rad2rad Valine{} = VAL+rad2rad Tryptophan{} = TRP+rad2rad Tyrosine{} = TYR
+ src/Bio/Protein/Chain.hs view
@@ -0,0 +1,25 @@+{-# OPTIONS_GHC -fno-warn-orphans #-}++module Bio.Protein.Chain+ ( module C+ , ProteinChain(..)+ ) where++import Bio.Chain as C+import Bio.Protein.AminoAcid+import Control.Lens+import Data.Array (Ix (..))+import Data.String (IsString (..))++newtype ProteinChain i a = ProteinChain { getChain :: Chain i a }+ deriving (Show, Eq, Functor, Foldable, Traversable, ChainLike)++type instance Index (ProteinChain i a) = i+type instance IxValue (ProteinChain i a) = a++instance Ix i => Ixed (ProteinChain i a) where+ ix i' = coerced . ix @(Chain i a) i'++instance IsString (ProteinChain Int AA) where+ fromString = ProteinChain . fromString+
+ src/Bio/Protein/Chain/Builder.hs view
@@ -0,0 +1,141 @@+{-# LANGUAGE TypeSynonymInstances #-}++module Bio.Protein.Chain.Builder+ ( Buildable (..)+ , build+ ) where++import Data.Ix ( Ix )+import Control.Lens+import Linear.V3 ( V3 (..)+ , cross+ , _z+ )+import Linear.Vector ( negated+ , unit+ , (*^)+ )++import Bio.Utils.Geometry hiding ( angle )+import Bio.Protein.AminoAcid+import Bio.Protein.Chain++class Buildable a where+ type Monomer a :: *+ initB :: Monomer a -> a+ nextB :: Monomer a -> a -> a++build :: forall a m.(Buildable a, ChainLike m, Ix (Index m), IxValue m ~ Monomer a) => m -> ProteinChain (Index m) a+build ch = ProteinChain result+ where+ result :: Chain (Index m) a+ result = chain (bounds ch) [ (i, next i x) | (i, x) <- assocs ch ]+ next :: Index m -> Monomer a -> a+ next k x | k == fst (bounds ch) = initB x+ | otherwise = nextB x (result ! pred k)++instance Buildable (BB V3R) where+ type Monomer (BB V3R) = AA++ -- | Place first amino acid backbone in some chain+ -- The placement will be like this:+ -- y /|\+ -- |+ -- |+ -- N | Ca+ -- ----*-----*------------->+ -- | C x+ -- | *+ -- |+ --+ initB _ = let n_ = V3 n_x 0.0 0.0+ a_ = V3 0.0 0.0 0.0+ c_ = V3 c_x c_y 0.0+ --+ n_x = - dist N CA+ c_x = dist CA C * cos (pi + angle N CA C)+ c_y = dist CA C * sin (pi + angle N CA C)+ in create @(BB V3R) n_ a_ c_++ -- | Place next amino acid backbone in some chain+ -- The placement can be done by two cases.+ -- First:+ -- Ca_i N_i+1 C_i+1+ -- * * *+ -- + -- * * *+ -- N_i C_i Ca_i+1+ -- Second:+ -- N_i C_i Ca_i+1+ -- * * *+ --+ -- * * *+ -- Ca_i N_i+1 C_i+1+ --+ -- Let us enumerate atoms: 1 for N_i, 2 for Ca_i, 3 for C_i, 4 for N_i+1, 5 for Ca_i+1, 6 for C_i+1.+ -- We have to find points 4, 5, 6 using 1, 2, 3. To find this points let us introduce vectors named+ -- like 'vij' from i to j, e.g. v12 is a vector from N_i to Ca_i. Our main idea will be to get a + -- direction vector from i+1 to i, rotate it and then upscale by specified bond length. One thing to+ -- look at is the direction of rotations. If we have the first case, then the first rotation should be+ -- conterclock-wise, otherwise — clock-wise. To detect it we have to understand whether 3 is on the left+ -- of 12 vector (first case) or on the right. We can understand it using v21 and v23:+ -- if (v21 `cross` v23) ^. _z < 0 then First else Second. First means that every angle should be negated.+ -- So, we can determine coordinate of 4. First we get the v32 and normalize it, then we will rotate it to+ -- CA-C-N angle (multiplied by -1 or not), next multiply this direction vector by typical C-N bond length+ -- and at last add the obtained vector to 3. The same idea is used to find point 5, but now we should+ -- make out rotation in the opposite direction. At last we will do the same with point 6.+ --+ nextB _ aa = let -- we will always rotate around Z+ rot = rotate (unit _z)+ -- determine the direction+ v21 = aa ^. n . atom - aa ^. ca . atom+ v23 = aa ^. c . atom - aa ^. ca . atom+ cw = if (v21 `cross` v23) ^. _z < 0 then 1.0 else -1.0 :: R+ -- determine the coordinate of n (point 4)+ v32 = negated v23+ v34 = dist C N *^ rot (cw * angle CA C N) (normalize v32)+ n_ = aa ^. c . atom + v34+ -- determine the coordinate of ca (point 5)+ v43 = negated v34+ v45 = dist N CA *^ rot (-cw * angle C N CA) (normalize v43)+ ca_ = n_ + v45+ -- determine the coordinate of ca (point 6)+ v54 = negated v45+ v56 = dist CA C *^ rot (cw * angle N CA C) (normalize v54)+ c_ = ca_ + v56+ in create @(BB V3R) n_ ca_ c_++instance Buildable (BBT V3R) where+ type Monomer (BBT V3R) = AA++ initB t = let aa = initB t :: BB V3R+ in create @(BBT V3R) (aa ^. n . atom) (aa ^. ca . atom) (aa ^. c . atom) t++ nextB t aaT = let aa = create @(BB V3R) (aaT ^. n . atom) (aaT ^. ca . atom) (aaT ^. c . atom)+ ab = nextB t aa :: BB V3R+ in create @(BBT V3R) (ab ^. n . atom) (ab ^. ca . atom) (ab ^. c . atom) t++-- Helper types and functions++-- | Atoms of amino acid backbone+--+data BackboneAtom = N | CA | C+ deriving (Show, Eq, Ord, Bounded, Enum)++-- | Atoms of amino acid radicals (TODO: fill this)+--+-- data RadicalAtom++-- | Distance between two basic backbone atom types+dist :: BackboneAtom -> BackboneAtom -> R+dist N CA = 1.460+dist CA C = 1.509+dist C N = 1.290+dist x y = dist y x++-- | Angles between every triple of succesive atoms+angle :: BackboneAtom -> BackboneAtom -> BackboneAtom -> R+angle N CA C = pi * 110.990 / 180.0+angle CA C N = pi * 118.995 / 180.0+angle C N CA = angle CA C N+angle x y z = angle z y x
+ src/Bio/Protein/Metric.hs view
@@ -0,0 +1,28 @@+module Bio.Protein.Metric+ ( Metricable (..)+ ) where++import Data.Monoid ( First (..) )+import Control.Lens+import Bio.Utils.Geometry ( V3R+ , R+ )+import qualified Bio.Utils.Geometry as G++class Metricable m where+ type ReturnMetric m :: *+ distance :: Getting m a V3R -> Getting m a V3R -> Getting (ReturnMetric m) a R+ angle :: Getting m a V3R -> Getting m a V3R -> Getting m a V3R -> Getting (ReturnMetric m) a R+ dihedral :: Getting m a V3R -> Getting m a V3R -> Getting m a V3R -> Getting m a V3R -> Getting (ReturnMetric m) a R++instance Metricable (First V3R) where+ type ReturnMetric (First V3R) = First R+ distance x y _ aa = Const . First $ G.distance <$> (aa ^? x) <*> (aa ^? y)+ angle x y z _ aa = Const . First $ G.angle <$> ((-) <$> aa ^? x <*> aa ^? y) <*> ((-) <$> aa ^? z <*> aa ^? y)+ dihedral x y z w _ aa = Const . First $ G.dihedral <$> (aa ^? x) <*> (aa ^? y) <*> (aa ^? z) <*> (aa ^? w)++instance Metricable V3R where+ type ReturnMetric V3R = R+ distance x y _ aa = Const $ G.distance (aa ^. x) (aa ^. y)+ angle x y z _ aa = Const $ G.angle (aa ^. x - aa ^. y) (aa ^. z - aa ^. y)+ dihedral x y z w _ aa = Const $ G.dihedral (aa ^. x) (aa ^. y) (aa ^. z) (aa ^. w)
+ src/Bio/Utils/Geometry.hs view
@@ -0,0 +1,108 @@+{-# LANGUAGE TemplateHaskell #-}++module Bio.Utils.Geometry+ ( R+ , V3R+ , Ray (..)+ , AffineTransformable(..)+ , Epsilon (..)+ , zoRay+ , cross, dot+ , norm , normalize+ , distance, angle, dihedral+ , svd3+ ) where++import Control.Lens+import Linear.V3 ( V3+ , cross+ )+import Linear.Vector ( zero )+import Linear.Epsilon ( Epsilon (..) )+import Linear.Matrix ( M33 )+import Linear.Metric ( dot+ , norm+ , normalize+ , distance+ )+import qualified Linear.Quaternion as Q+ ( rotate+ , axisAngle+ )++-- | Default floating point type, switch here to move to Doubles+--+type R = Float++-- | Defalut type of 3D vectors+--+type V3R = V3 R++-- | Ray has an origin and a direction+--+data Ray a = Ray { _origin :: a+ , _direction :: a+ }++makeLenses ''Ray++-- | Zero-origin ray+zoRay :: V3R -> Ray V3R+zoRay = Ray zero . normalize++-- | Affine transformations for vectors and sets of vectors+--+class AffineTransformable a where+ -- | Rotate an object around the vector by some angle+ --+ rotate :: V3R -> R -> a -> a++ -- | Rotate an object around the ray by some angle+ --+ rotateR :: Ray V3R -> R -> a -> a++ -- | Translocate an object by some vectors+ --+ translate :: V3R -> a -> a++-- | We can apply affine transformations to vectors+--+instance AffineTransformable V3R where+ rotate v a = Q.rotate (Q.axisAngle v a)+ rotateR r a x = rotate (r ^. direction) a (x - r ^. origin) + r ^. origin+ translate v = (v +)++-- | If we have any collection of vectors, than we can transform it too+--+instance Functor f => AffineTransformable (f V3R) where+ rotate v a = fmap (rotate v a)+ rotateR r a = fmap (rotateR r a)+ translate v = fmap (translate v)++-- | Measure angle between vectors+--+angle :: V3R -> V3R -> R+angle a b = atan2 (norm (a `cross` b)) (a `dot` b)++-- | Measure dihedral between four points+-- by https://math.stackexchange.com/a/47084+--+dihedral :: V3R -> V3R -> V3R -> V3R -> R+dihedral x y z w = let b1 = y - x+ b2 = z - y+ b3 = w - z+ n1 = normalize $ b1 `cross` b2+ n2 = normalize $ b2 `cross` b3+ m1 = n1 `cross` normalize b2+ in atan2 (m1 `dot` n2) (n1 `dot` n2)++data SVD a = SVD { svdU :: a+ , svdS :: a+ , svdV :: a+ }+ deriving (Show, Eq)++-- | Singular value decomposition+-- for 3x3 matricies+svd3 :: M33 R -> SVD (M33 R)+svd3 = undefined
+ src/Bio/Utils/IUPAC.hs view
@@ -0,0 +1,39 @@+module Bio.Utils.IUPAC+ (+ AtomType (..)+ ) where++import Control.DeepSeq (NFData (..))+import GHC.Generics (Generic)++-- | Atom types in IUPAC nomenclature.+--+data AtomType = N | CA | C | O | OXT+ | CB+ | CG | CG1 | CG2+ | CD | CD1 | CD2+ | CE | CE1 | CE2 | CE3+ | CH3+ | CZ | CZ2 | CZ3+ | CH2+ | SG+ | SD+ | OG | OG1+ | OD1 | OD2+ | OE1 | OE2+ | OH+ | ND1 | ND2+ | NE | NE1 | NE2+ | NZ+ | NH1 | NH2+ | H+ | HA | HA2 | HA3+ | HB | HB1 | HB2 | HB3+ | HG | HG1 | HG2 | HG3 | HG11 | HG12 | HG13 | HG21 | HG22 | HG23+ | HD | HD1 | HD2 | HD3 | HD11 | HD12 | HD13 | HD21 | HD22 | HD23+ | HE | HE1 | HE2 | HE3 | HE21 | HE22+ | HH | HH2 | HH11 | HH12 | HH21 | HH22+ | HZ | HZ1 | HZ2 | HZ3+ deriving (Show, Read, Eq, Ord, Generic)++instance NFData AtomType
+ src/Bio/Utils/Matrix.hs view
@@ -0,0 +1,20 @@+module Bio.Utils.Matrix+ ( eig3+ ) where++import Linear.Matrix+import Linear.V3++import Bio.Utils.Geometry ( R )++eig3 :: M33 R -> V3 R+eig3 m | isSym m = ei3sym m+ | otherwise = undefined+ where+ ei3sym :: M33 R -> V3 R+ ei3sym = undefined++isSym :: Eq a => M33 a -> Bool+isSym (V3 (V3 _ a21 a31)+ (V3 a12 _ a32)+ (V3 a13 a23 _ )) = a21 == a12 && a31 == a13 && a32 == a23
+ src/Bio/Utils/Monomer.hs view
@@ -0,0 +1,30 @@+module Bio.Utils.Monomer+ ( Symbol(..)+ , FromSymbol(..)+ , ThreeSymbols(..)+ , FromThreeSymbols (..)+ ) where++import Data.Text ( Text )++class Symbol a where+ symbol :: a -> Char++instance Symbol Char where+ symbol = id++class FromSymbol a where+ fromSymbol :: Char -> Maybe a+ fromSymbol = either (const Nothing) Just . fromSymbolE++ fromSymbolE :: Char -> Either Char a+ fromSymbolE c = maybe (Left c) Right $ fromSymbol c++instance FromSymbol Char where+ fromSymbolE = Right++class ThreeSymbols a where+ threeSymbols :: a -> Text++class FromThreeSymbols a where+ fromThreeSymbols :: Text -> Maybe a
+ test/HandcraftedSpec.hs view
@@ -0,0 +1,57 @@+module HandcraftedSpec where++import Bio.Chain.Alignment+import Bio.Chain.Alignment.Scoring+import Test.Hspec++handcraftedTests :: Spec+handcraftedTests = do+ handcraftedLocal++handcraftedLocal :: Spec+handcraftedLocal = do+ describe "LocalAlignment. AffineGap2" $ do+ let testAlign = align (LocalAlignment nuc44 (AffineGap (-5) (-1), AffineGap (-1000) (-1000)))+ let a = "AAAAAAGGGGGGGGGGGGTTTTTTTTT"+ b = "AAAAAATTTTTTTTT"+ aAns = "TTTTTTTTT"+ bAns = "TTTTTTTTT"+ res = testAlign (a :: String) (b :: String)+ (a', b') = viewAlignment res+ score1 = score res+ it "first sequence" $ a' `shouldBe` aAns+ it "second sequence" $ b' `shouldBe` bAns+ it "score" $ score1 `shouldBe` 45++ let a2 = "AAAAAATTTTTTTTT"+ b2 = "AAAAAAGGGGGGGGGGGGTTTTTTTTT"+ a2Ans = "AAAAAA------------TTTTTTTTT"+ b2Ans = "AAAAAAGGGGGGGGGGGGTTTTTTTTT"+ res2 = testAlign (a2 :: String) (b2 :: String)+ (a2', b2') = viewAlignment res2+ score2 = score res2+ it "first sequence" $ a2' `shouldBe` a2Ans+ it "second sequence" $ b2' `shouldBe` b2Ans+ it "score" $ score2 `shouldBe` 59++ let a3 = "AACCTTCCAACCGGCCAACCTTCCAACCGGCCTT"+ b3 = "AATTAAGGAATTAAGGTT"+ a3Ans = "GGCCTT"+ b3Ans = "GGAATT"+ res3 = testAlign (a3 :: String) (b3 :: String)+ (a3', b3') = viewAlignment res3+ score3 = score res3+ it "first sequence" $ a3' `shouldBe` a3Ans+ it "second sequence" $ b3' `shouldBe` b3Ans+ it "score" $ score3 `shouldBe` 12++ let a4 = "AATTAAGGAATTAAGGTT"+ b4 = "AACCTTCCAACCGGCCAACCTTCCAACCGGCCTT"+ a4Ans = "AA--TT--AA--GG--AA--TT--AA--GG--TT"+ b4Ans = "AACCTTCCAACCGGCCAACCTTCCAACCGGCCTT"+ res4 = testAlign (a4 :: String) (b4 :: String)+ (a4', b4') = viewAlignment res4+ score4 = score res4+ it "first sequence" $ a4' `shouldBe` a4Ans+ it "second sequence" $ b4' `shouldBe` b4Ans+ it "score" $ score4 `shouldBe` 42
+ test/JuliaSpec.hs view
@@ -0,0 +1,270 @@+module JuliaSpec where++import Bio.Chain.Alignment+import Bio.Chain.Alignment.Scoring+import Test.Hspec++juliaBasedTests :: Spec+juliaBasedTests = do+ juliaGlobal+ juliaSemiglobal+ juliaLocal++juliaGlobal :: Spec+juliaGlobal = do+ describe "Global Alignment, Simple Gap" $ do+ let testAlign = align (GlobalAlignment nuc44 (-9 :: Int))++ let resBig1 = testAlign (bigA :: String) (bigB :: String)+ resBig2 = testAlign (bigB :: String) (bigA :: String)+ resMed1 = testAlign (bigA :: String) (bigB :: String)+ resMed2 = testAlign (bigB :: String) (bigA :: String)+ resSmall1 = testAlign (bigA :: String) (bigB :: String)+ resSmall2 = testAlign (bigB :: String) (bigA :: String)++ it "Is symmetric (big): " $ score resBig1 `shouldBe` score resBig2+ it "Is symmetric (med): " $ score resMed1 `shouldBe` score resMed2+ it "Is symmetric (small): " $ score resSmall1 `shouldBe` score resSmall2+ -- | scoremodel = AffineGapScoreModel(EDNAFULL, gap_open=0, gap_extend=-9);+ -- | pairalign(GlobalAlignment(), a, b, scoremodel)+ let aAns = "CTCTGCAAGCATAGAGATATTTCACCGGCAATATTTTTCGTTGAAGTGTATTTGTCCCTATTATCACACC\+ \AGTCATTCCGATGTCGTTAGGGCCGCCTGTTTTACTGGAGTTCCCGTTCGATGTACCCTCCTCCATAGA\+ \GATTTAGGGGTATAGACCGCGACGGACGGGGCTTGCGAAAGCTGCCCGATAATGGACTTTGACAATAGA\+ \CGTAATCCGTGAACGGCTGGCGTTTCACACTCAACCGCATTGAACAGCCATGCCATCCGCAGGATGCGC\+ \CAGGGAAAAGCCACCCCCAAGAACGACTCCAGGGGCCCAAGTATGATGATAGCGTTAGGCTTGCTCGAA\+ \GGCAAACAACTCCGAGGATGCTTAGCAGGCCACAACGTGCTCGATAGCACGAGCGTAAGATGCGAACGC\+ \AGTAGTACCTAATGCGGGCGAGACATTTCGGGCCGTTGAAGCCCTGTCTGTTTTTTCGTACAGAGGGTGGTTTGAAGTATCGCCC"+ bAns = "-T--G-AAGCA--GC-AT-TCG-ACT--CACTCATGT-CGCCTA-GAG-A---GAA--T-TTAT-A----\+ \-GTCAA----ATGAC-TT-GTG--GC--GTTCT--TGGACTT----TTAGA-G--CCGAAATGCAT-G-\+ \GA----GG---A-A-AC----A----------TT---AA-G--G----ATG-T--A-TTTG----T---\+ \--TA--C--TGA--G-C---------A---T-----G-ATT-AA--GC-A----AT--------TGC-C\+ \-AGGTAAA--C----------------T---G---CC-A-G---G-T--T----T--GG--------A-\+ \--C---CAA-TC------T--T--G--GG--A----G--CT---TAGTAC---C-TAA-A-GA-AA-G-\+ \A--AGT-C--A-TGCA--C-AGTCA---CG--CCGC--A--CC----------------A-A-----T------A--TA-CAC--"+ res = testAlign (bigA :: String) (bigB :: String)+ (a', b') = viewAlignment res+ score1 = score res+ it "first sequence" $ a' `shouldBe` aAns+ it "second sequence" $ b' `shouldBe` bAns+ it "score" $ score1 `shouldBe` (-1605)++ let a2Ans = "-GT-TTGGG-CATATTC-AAG-ATCA-GAC-C-A----ATCGGTC-GAT-G--TGAAACAAGT--TA-\+ \---ATAATGCTA-CACAGTGTTC--GCTG-T-TTT-A-CTTCGG-GT-CCTCCTGCCCCTTGAGG-A\+ \--TATGAT--A--CAC----AACCG--CT--CTCTAGACA-GGAAGAAA-ACTGCACCGGATA--"+ b2Ans = "AGCATTAGCTCACATAATAAGGATTGTGGGGCTACGGCATTTGTATGATAGCCTGACCGAGGTGGTAG\+ \GGCATTAAGGTATCGTCGTCCTCCAGCTGGTGTTTTATCTGCCGCGCACCTCAAGTCGACTACGGTA\+ \ACTATGGTTTATGCACGTGCAAACGTACTAACACTGGAGGCGGAATAATGATTAGAGTGGCCCGT"+ res2 = testAlign (medA :: String) (medB :: String)+ (a2', b2') = viewAlignment res2+ score2 = score res2++ it "first sequence" $ a2' `shouldBe` a2Ans+ it "second sequence" $ b2' `shouldBe` b2Ans+ it "score" $ score2 `shouldBe` (-168)++ let a3Ans = "ATCATCTCGGACTATGGTTGGGGATCTAAGACGTCCTCTTGGCTCTCGCC"+ b3Ans = "AAGAGATCGT--TATCGC-GCGGAT-TAAT---TCCGAGTCG-TC-CAC-"+ res3 = testAlign (smallA :: String) (smallB :: String)+ (a3', b3') = viewAlignment res3+ score3 = score res3++ it "first sequence" $ a3' `shouldBe` a3Ans+ it "second sequence" $ b3' `shouldBe` b3Ans+ it "score" $ score3 `shouldBe` (-16)++ describe "GlobalAlignment. AffineGap" $ do+ let testAlign = align (GlobalAlignment nuc44 (AffineGap (-5) (-1)))+ res1 = testAlign "AGGT" "CAAT"+ res2 = testAlign "AGGTACC" "CAATG"++ -- Честно посчитано руками в тетради+ it "score" $ score res1 `shouldBe` (-4)+ it "score" $ score res2 `shouldBe` (-8)+ it "alignment" $ viewAlignment res1 `shouldBe` ("-AGGT", "CAA-T")+ it "alignment" $ viewAlignment res2 `shouldBe` ("--AGGTACC","CAATG----")++ describe "GlobalAlignment. Sanity test" $ do+ let resBig1 = align (GlobalAlignment nuc44 (AffineGap (-5) (-5))) (bigA :: String) (bigB :: String)+ resBig2 = align (GlobalAlignment nuc44 (-5 :: Int)) (bigB :: String) (bigA :: String)++ resMed1 = align (GlobalAlignment nuc44 (AffineGap (-5) (-5))) (medA :: String) (medB :: String)+ resMed2 = align (GlobalAlignment nuc44 (-5 :: Int)) (medB :: String) (medA :: String)++ resSmall1 = align (GlobalAlignment nuc44 (AffineGap (-5) (-5))) (smallA :: String) (smallB :: String)+ resSmall2 = align (GlobalAlignment nuc44 (-5 :: Int)) (smallB :: String) (smallA :: String)+ it "SimpleGap is a special case of AffineGap (big)" $ score resBig1 `shouldBe` score resBig2+ it "SimpleGap is a special case of AffineGap (med)" $ score resMed1 `shouldBe` score resMed2+ it "SimpleGap is a special case of AffineGap (small)" $ score resSmall1 `shouldBe` score resSmall2+++juliaSemiglobal :: Spec+juliaSemiglobal = do+ describe "Semiglobal Alignment, Simple Gap" $ do+ let testAlign = align (SemiglobalAlignment nuc44 (-3 :: Int))+ let resBig1 = testAlign (bigA :: String) (bigB :: String)+ resBig2 = testAlign (bigB :: String) (bigA :: String)+ resMed1 = testAlign (bigA :: String) (bigB :: String)+ resMed2 = testAlign (bigB :: String) (bigA :: String)+ resSmall1 = testAlign (bigA :: String) (bigB :: String)+ resSmall2 = testAlign (bigB :: String) (bigA :: String)++ it "Is symmetric (big): " $ score resBig1 `shouldBe` score resBig2+ it "Is symmetric (med): " $ score resMed1 `shouldBe` score resMed2+ it "Is symmetric (small): " $ score resSmall1 `shouldBe` score resSmall2+ -- | scoremodel = AffineGapScoreModel(EDNAFULL, gap_open=0, gap_extend=-3);+ -- | pairalign(OverlapAlignment(), a, b, scoremodel)+ let aAns = "CTCTGCAAGCATAGAGATATTTCACCGGCAATATTTTTCGTTGAAGTGTATTTGTCCCTATTATCACACC\+ \AGTCATTCCGATGTCGTTAGGGCCGCCTGTTTTACTG--G-AG--TTC--C-CGTTCGATGTACCCTCC\+ \TCCATAGAGA-TTTAGGGGTA-T-AGACCGCGACGGACGGGGCTTGCGAAAGCTGCCCGATAATGGACT\+ \TTGACAATAGACGTAATCCGTGAACGGC-TGGCGTTTCACACTCAACCGCATTGAACAGCCATGCCATC\+ \CGC-AGGATGCGCCAGGGAA-AAGCCACCCCCAAG-AACGACT-CCAGG---GGCCCAAGTATGATGAT\+ \AGCGTTAGGCTTGCTCGAAGGCAAACAACTCCGAGGATGCTTAG-CAGGCCACAACGTGC--TCGATAG\+ \CACGAGCGTAAGATGCGAACGCAGTAGTACCTAATGCGGGCGAGACATTTCGGGCCGTTGAAGCCCTGT\+ \CTGTTTTTTCGTACAGAGGGTGGTTTGAAGTATCGCCC"+ bAns = "----------------------------------------------------------------------\+ \-----------------------------------TGAAGCAGCATTCGACTCACTC-ATGT-CGC-C-\+ \T--AGAGAGAATTTA----TAGTCA-A-----ATG-AC-----TTGTG---GC-GTTC--T--TGGACT\+ \TT-----TAGA-G----CCGA-AATG-CATGGAGG---A-A---A-C---ATT-AAG-G--ATGT-ATT\+ \TGTTAC--TGAGC-ATG-ATTAAGCAATTGCCAGGTAA--ACTGCCAGGTTTGGACCAA-TCT--TGGG\+ \AGC-TTAG--TACCTA-AAG--AAAGAAGTC--A---TGCACAGTCACGCCGCA-C---CAAT--ATA-\+ \CAC------------------------------------------------------------------\+ \--------------------------------------"+ res = testAlign (bigA :: String) (bigB :: String)+ (a', b') = viewAlignment res+ score1 = score res+ it "first sequence" $ a' `shouldBe` aAns+ it "second sequence" $ b' `shouldBe` bAns+ it "score" $ score1 `shouldBe` 358++ let a2Ans = "---------------------G-TT-TGGG-C-AT---ATTCA-A-GATCAGACCAATCGGTCGATGT\+ \GAAA---CAAGTTAA--TAAT-G-C-TACAC-AG-TGTTCGCTGTTTTACTTCGGGTCCTCCTGCCC\+ \CTTGAG--GA-TATGATACACAACCGCTCTCTA-G-ACAG-G-AA--G-A--AA-ACTGCAC-CGGA\+ \-TA------------------"+ b2Ans = "AGCATTAGCTCACATAATAAGGATTGTGGGGCTACGGCATTTGTATGAT-AG-CC--T-GACCGAGGT\+ \GGTAGGGCA--TTAAGGTA-TCGTCGTCCTCCAGCTG---G-TGTTTTA-T-CTG--CCGC--GCAC\+ \CTCAAGTCGACTACGGTA-ACTATGG-T-T-TATGCAC-GTGCAAACGTACTAACACTGGAGGCGGA\+ \ATAATGATTAGAGTGGCCCGT"+ res2 = testAlign (medA :: String) (medB :: String)+ (a2', b2') = viewAlignment res2+ score2 = score res2++ it "first sequence" $ a2' `shouldBe` a2Ans+ it "second sequence" $ b2' `shouldBe` b2Ans+ it "score" $ score2 `shouldBe` 266++ let a3Ans = "ATCATCTCG--GACTATGGTT---G-G-GGATCTAA----GA--CGTCCTCTTGGCTCTCGCC"+ b3Ans = "---------AAGAG-ATCGTTATCGCGCGGAT-TAATTCCGAGTCGTCCAC------------"+ res3 = testAlign (smallA :: String) (smallB :: String)+ (a3', b3') = viewAlignment res3+ score3 = score res3++ it "first sequence" $ a3' `shouldBe` a3Ans+ it "second sequence" $ b3' `shouldBe` b3Ans+ it "score" $ score3 `shouldBe` 63+ describe "SemiglobalAlignment. Sanity test" $ do+ let resBig1 = align (SemiglobalAlignment nuc44 (AffineGap (-3) (-3))) (bigA :: String) (bigB :: String)+ resBig2 = align (SemiglobalAlignment nuc44 (-3 :: Int)) (bigB :: String) (bigA :: String)++ resMed1 = align (SemiglobalAlignment nuc44 (AffineGap (-3) (-3))) (medA :: String) (medB :: String)+ resMed2 = align (SemiglobalAlignment nuc44 (-3 :: Int)) (medB :: String) (medA :: String)++ resSmall1 = align (SemiglobalAlignment nuc44 (AffineGap (-3) (-3))) (smallA :: String) (smallB :: String)+ resSmall2 = align (SemiglobalAlignment nuc44 (-3 :: Int)) (smallB :: String) (smallA :: String)+ it "SimpleGap is a special case of AffineGap (big)" $ score resBig1 `shouldBe` score resBig2+ it "SimpleGap is a special case of AffineGap (med)" $ score resMed1 `shouldBe` score resMed2+ it "SimpleGap is a special case of AffineGap (small)" $ score resSmall1 `shouldBe` score resSmall2++juliaLocal :: Spec+juliaLocal = do+ describe "Local Alignment, Simple Gap" $ do+ let testAlign = align (LocalAlignment nuc44 (-7 :: Int))+ let resBig1 = testAlign (bigA :: String) (bigB :: String)+ resBig2 = testAlign (bigB :: String) (bigA :: String)+ resMed1 = testAlign (bigA :: String) (bigB :: String)+ resMed2 = testAlign (bigB :: String) (bigA :: String)+ resSmall1 = testAlign (bigA :: String) (bigB :: String)+ resSmall2 = testAlign (bigB :: String) (bigA :: String)++ it "Is symmetric (big): " $ score resBig1 `shouldBe` score resBig2+ it "Is symmetric (med): " $ score resMed1 `shouldBe` score resMed2+ it "Is symmetric (small): " $ score resSmall1 `shouldBe` score resSmall2+ -- | scoremodel = AffineGapScoreModel(EDNAFULL, gap_open=0, gap_extend=-7);+ -- | pairalign(LocalAlignment(), a, b, scoremodel)+ let aAns = "TGAA-CAGCCATGCCA-TC-CGCAGGATGCGCC-AGGGAAAAGCCACCCCCAAGAACGACTCCA\+ \GGGGCCCAAGTATGATGATAGC-GTTAGGCTTGCTCGAAGGCAAACAACTCCGAGG-ATGC-T\+ \TAGCAGGCCACAACGTGCTCGATAGCACGAGCGTAAGATGCGAACGCAGTAGTACCTAATGCGGGCGAGACATT"+ bAns = "TGAAGCAGC-ATTCGACTCACTCATG-T-CGCCTAGAGAGAATTTATAGTCAA-AT-GACTTGT\+ \GGCGTTCTTGGACTTTTAGAGCCGAAATGCATG---GA-GG-AAACATTAAGGATGTATTTGT\+ \TA-CTGAGCATGAT-TAAGCAATTGC-C-AG-GTAAACTGCCAG-GTT-TGG-ACC-AATCTTGG-GAG-C-TT"+ res = testAlign (bigA :: String) (bigB :: String)+ (a', b') = viewAlignment res+ score1 = score res+ it "first sequence" $ a' `shouldBe` aAns+ it "second sequence" $ b' `shouldBe` bAns+ it "score" $ score1 `shouldBe` 108++ let a2Ans = "GGCATATTCAAGATCAGACCAATCG-GTCGAT-GTGAAACAAGTTAATAATGCTACACAGTGT\+ \TCGCTGTTTTA-CTTCGG-GT-CCTCCTGCCCCTTGAGG-A--TATGATACACAACCGCTCT\+ \CTAGACAGGAAGAAAACTGCAC-CGGA-TA"+ b2Ans = "GGCATTTGTATGAT-AGCCTGACCGAGGTGGTAGGGCATTAAGGTA-TCGTCGTCCTCCA-GC\+ \TGG-TGTTTTATCTGCCGCGCACCTCAAGTCGACTACGGTAACTATGGTTTATGCACG-TG-\+ \CAA-AC-GTACTAACACTGGAGGCGGAATA"+ res2 = testAlign (medA :: String) (medB :: String)+ (a2', b2') = viewAlignment res2+ score2 = score res2++ it "first sequence" $ a2' `shouldBe` a2Ans+ it "second sequence" $ b2' `shouldBe` b2Ans+ it "score" $ score2 `shouldBe` 91++ let a3Ans = "TCTAAGACGTCCTC"+ b3Ans = "TCCGAGTCGTCCAC"+ res3 = testAlign (smallA :: String) (smallB :: String)+ (a3', b3') = viewAlignment res3+ score3 = score res3++ it "first sequence" $ a3' `shouldBe` a3Ans+ it "second sequence" $ b3' `shouldBe` b3Ans+ it "score" $ score3 `shouldBe` 34++ describe "LocalAlignment. Sanity test" $ do+ let resBig1 = align (LocalAlignment nuc44 (AffineGap (-7) (-7))) (bigA :: String) (bigB :: String)+ resBig2 = align (LocalAlignment nuc44 (-7 :: Int)) (bigB :: String) (bigA :: String)++ resMed1 = align (LocalAlignment nuc44 (AffineGap (-7) (-7))) (medA :: String) (medB :: String)+ resMed2 = align (LocalAlignment nuc44 (-7 :: Int)) (medB :: String) (medA :: String)++ resSmall1 = align (LocalAlignment nuc44 (AffineGap (-7) (-7))) (smallA :: String) (smallB :: String)+ resSmall2 = align (LocalAlignment nuc44 (-7 :: Int)) (smallB :: String) (smallA :: String)+ it "SimpleGap is a special case of AffineGap (big)" $ score resBig1 `shouldBe` score resBig2+ it "SimpleGap is a special case of AffineGap (med)" $ score resMed1 `shouldBe` score resMed2+ it "SimpleGap is a special case of AffineGap (small)" $ score resSmall1 `shouldBe` score resSmall2++bigA :: String+bigA = "CTCTGCAAGCATAGAGATATTTCACCGGCAATATTTTTCGTTGAAGTGTATTTGTCCCTATTATCACACCAG\+ \TCATTCCGATGTCGTTAGGGCCGCCTGTTTTACTGGAGTTCCCGTTCGATGTACCCTCCTCCATAGAGATT\+ \TAGGGGTATAGACCGCGACGGACGGGGCTTGCGAAAGCTGCCCGATAATGGACTTTGACAATAGACGTAAT\+ \CCGTGAACGGCTGGCGTTTCACACTCAACCGCATTGAACAGCCATGCCATCCGCAGGATGCGCCAGGGAAA\+ \AGCCACCCCCAAGAACGACTCCAGGGGCCCAAGTATGATGATAGCGTTAGGCTTGCTCGAAGGCAAACAAC\+ \TCCGAGGATGCTTAGCAGGCCACAACGTGCTCGATAGCACGAGCGTAAGATGCGAACGCAGTAGTACCTAA\+ \TGCGGGCGAGACATTTCGGGCCGTTGAAGCCCTGTCTGTTTTTTCGTACAGAGGGTGGTTTGAAGTATCGCCC"++bigB :: String+bigB = "TGAAGCAGCATTCGACTCACTCATGTCGCCTAGAGAGAATTTATAGTCAAATGACTTGTGGCGTTCTTGGACT\+ \TTTAGAGCCGAAATGCATGGAGGAAACATTAAGGATGTATTTGTTACTGAGCATGATTAAGCAATTGCCAGG\+ \TAAACTGCCAGGTTTGGACCAATCTTGGGAGCTTAGTACCTAAAGAAAGAAGTCATGCACAGTCACGCCGCACCAATATACAC"++medA :: String+medA = "GTTTGGGCATATTCAAGATCAGACCAATCGGTCGATGTGAAACAAGTTAATAATGCTACACAGTGTTCGCTG\+ \TTTTACTTCGGGTCCTCCTGCCCCTTGAGGATATGATACACAACCGCTCTCTAGACAGGAAGAAAACTGCACCGGATA"++medB :: String+medB = "AGCATTAGCTCACATAATAAGGATTGTGGGGCTACGGCATTTGTATGATAGCCTGACCGAGGTGGTAGGGCA\+ \TTAAGGTATCGTCGTCCTCCAGCTGGTGTTTTATCTGCCGCGCACCTCAAGTCGACTACGGTAACTATGGT\+ \TTATGCACGTGCAAACGTACTAACACTGGAGGCGGAATAATGATTAGAGTGGCCCGT"++smallA :: String+smallA = "ATCATCTCGGACTATGGTTGGGGATCTAAGACGTCCTCTTGGCTCTCGCC"++smallB :: String+smallB = "AAGAGATCGTTATCGCGCGGATTAATTCCGAGTCGTCCAC"
+ test/Spec.hs view
@@ -0,0 +1,87 @@+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE TypeApplications #-}++import HandcraftedSpec (handcraftedTests)+import JuliaSpec (juliaBasedTests)+import Test.Hspec++import Control.Lens+import Data.Text++import Bio.Protein.AminoAcid+import Bio.Protein.Chain hiding (chain)+import Bio.Protein.Chain.Builder+import Bio.Utils.Geometry++buildChainSpec :: Spec+buildChainSpec = describe "Chain builder (BBT)" $ do+ let chain = build [ALA, CYS, ASP] :: ProteinChain Int (BBT V3R)+ let aa1 = chain ^?! ix 0+ aa2 = chain ^?! ix 1+ aa3 = chain ^?! ix 2+ let nac = pi * 110.990 / 180.0+ acn = pi * 118.995 / 180.0+ cna = pi * 118.995 / 180.0+ na = 1.460+ ac = 1.509+ cn = 1.290+ let a_ = aa1 ^. ca . atom+ n_ = aa1 ^. n . atom+ c_ = aa1 ^. c . atom+ let a2_ = aa2 ^. ca . atom+ n2_ = aa2 ^. n . atom+ c2_ = aa2 ^. c . atom+ let a3_ = aa3 ^. ca . atom+ n3_ = aa3 ^. n . atom+ c3_ = aa3 ^. c . atom+ it "builds single amino acid" $ do+ distance n_ a_ - na `shouldSatisfy` nearZero+ distance a_ c_ - ac `shouldSatisfy` nearZero+ angle (a_ - n_) (a_ - c_) - nac `shouldSatisfy` nearZero+ it "builds even amino acids" $ do+ distance n2_ c_ - cn `shouldSatisfy` nearZero+ angle (c_ - a_) (c_ - n2_) - acn `shouldSatisfy` nearZero+ angle (n2_ - c_) (n2_ - a2_) - cna `shouldSatisfy` nearZero+ distance n2_ a2_ - na `shouldSatisfy` nearZero+ distance a2_ c2_ - ac `shouldSatisfy` nearZero+ angle (a2_ - n2_) (a2_ - c2_) - nac `shouldSatisfy` nearZero+ it "builds odd amino acids" $ do+ distance n3_ c2_ - cn `shouldSatisfy` nearZero+ angle (c2_ - a2_) (c2_ - n3_) - acn `shouldSatisfy` nearZero+ angle (n3_ - c2_) (n3_ - a3_) - cna `shouldSatisfy` nearZero+ distance n3_ c2_ - cn `shouldSatisfy` nearZero+ distance n3_ a3_ - na `shouldSatisfy` nearZero+ distance a3_ c3_ - ac `shouldSatisfy` nearZero+ angle (a3_ - n3_) (a3_ - c3_) - nac `shouldSatisfy` nearZero++lensesSpec :: Spec+lensesSpec = describe "Amino acid lenses" $ do+ it "works on BB" $ do+ let aa = create @(BB Text) "N" "CA" "C"+ aa ^. n . atom `shouldBe` "N"+ aa ^. ca . atom `shouldBe` "CA"+ aa ^. c . atom `shouldBe` "C"+ it "works on BBT" $ do+ let aa = create @(BBT Text) "N" "CA" "C" ALA+ aa ^. n . atom `shouldBe` "N"+ aa ^. ca . atom `shouldBe` "CA"+ aa ^. c . atom `shouldBe` "C"+ aa ^. radical `shouldBe` ALA+ it "works on BBO" $ do+ let aa = create @(BBO Text) "N" "CA" "C" "O"+ aa ^. n . atom `shouldBe` "N"+ aa ^. ca . atom `shouldBe` "CA"+ aa ^. c . atom `shouldBe` "C"+ aa ^. o . atom `shouldBe` "O"+++alignmentSpec :: Spec+alignmentSpec = describe "Alignment" $ do+ juliaBasedTests+ handcraftedTests++main :: IO ()+main = hspec $ do+ lensesSpec+ buildChainSpec+ alignmentSpec