diff --git a/ChangeLog.md b/ChangeLog.md new file mode 100644 --- /dev/null +++ b/ChangeLog.md @@ -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. diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,4 @@ +# cobot + +[![Travis](https://img.shields.io/travis/less-wrong/cobot.svg)](https://travis-ci.org/less-wrong/cobot) +[![license](https://img.shields.io/github/license/less-wrong/cobot.svg)]() diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/bench/Main.hs b/bench/Main.hs new file mode 100644 --- /dev/null +++ b/bench/Main.hs @@ -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 + ] + ] diff --git a/cobot.cabal b/cobot.cabal new file mode 100644 --- /dev/null +++ b/cobot.cabal @@ -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 +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 diff --git a/src/Bio/Chain.hs b/src/Bio/Chain.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Chain.hs @@ -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 diff --git a/src/Bio/Chain/Alignment.hs b/src/Bio/Chain/Alignment.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Chain/Alignment.hs @@ -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 diff --git a/src/Bio/Chain/Alignment/Algorithms.hs b/src/Bio/Chain/Alignment/Algorithms.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Chain/Alignment/Algorithms.hs @@ -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 diff --git a/src/Bio/Chain/Alignment/Scoring.hs b/src/Bio/Chain/Alignment/Scoring.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Chain/Alignment/Scoring.hs @@ -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 +|] diff --git a/src/Bio/Chain/Alignment/Scoring/Loader.hs b/src/Bio/Chain/Alignment/Scoring/Loader.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Chain/Alignment/Scoring/Loader.hs @@ -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 diff --git a/src/Bio/Chain/Alignment/Scoring/TH.hs b/src/Bio/Chain/Alignment/Scoring/TH.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Chain/Alignment/Scoring/TH.hs @@ -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 diff --git a/src/Bio/Chain/Alignment/Type.hs b/src/Bio/Chain/Alignment/Type.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Chain/Alignment/Type.hs @@ -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' diff --git a/src/Bio/Molecule.hs b/src/Bio/Molecule.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Molecule.hs @@ -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 diff --git a/src/Bio/NucleicAcid/Chain.hs b/src/Bio/NucleicAcid/Chain.hs new file mode 100644 --- /dev/null +++ b/src/Bio/NucleicAcid/Chain.hs @@ -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 diff --git a/src/Bio/NucleicAcid/Nucleotide.hs b/src/Bio/NucleicAcid/Nucleotide.hs new file mode 100644 --- /dev/null +++ b/src/Bio/NucleicAcid/Nucleotide.hs @@ -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 diff --git a/src/Bio/NucleicAcid/Nucleotide/Instances.hs b/src/Bio/NucleicAcid/Nucleotide/Instances.hs new file mode 100644 --- /dev/null +++ b/src/Bio/NucleicAcid/Nucleotide/Instances.hs @@ -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 diff --git a/src/Bio/NucleicAcid/Nucleotide/Type.hs b/src/Bio/NucleicAcid/Nucleotide/Type.hs new file mode 100644 --- /dev/null +++ b/src/Bio/NucleicAcid/Nucleotide/Type.hs @@ -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 + diff --git a/src/Bio/Protein/Algebra.hs b/src/Bio/Protein/Algebra.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Protein/Algebra.hs @@ -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 diff --git a/src/Bio/Protein/AminoAcid.hs b/src/Bio/Protein/AminoAcid.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Protein/AminoAcid.hs @@ -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 diff --git a/src/Bio/Protein/AminoAcid/Instances.hs b/src/Bio/Protein/AminoAcid/Instances.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Protein/AminoAcid/Instances.hs @@ -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 diff --git a/src/Bio/Protein/AminoAcid/Type.hs b/src/Bio/Protein/AminoAcid/Type.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Protein/AminoAcid/Type.hs @@ -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 diff --git a/src/Bio/Protein/Chain.hs b/src/Bio/Protein/Chain.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Protein/Chain.hs @@ -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 + diff --git a/src/Bio/Protein/Chain/Builder.hs b/src/Bio/Protein/Chain/Builder.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Protein/Chain/Builder.hs @@ -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 diff --git a/src/Bio/Protein/Metric.hs b/src/Bio/Protein/Metric.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Protein/Metric.hs @@ -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) diff --git a/src/Bio/Utils/Geometry.hs b/src/Bio/Utils/Geometry.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Utils/Geometry.hs @@ -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 diff --git a/src/Bio/Utils/IUPAC.hs b/src/Bio/Utils/IUPAC.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Utils/IUPAC.hs @@ -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 diff --git a/src/Bio/Utils/Matrix.hs b/src/Bio/Utils/Matrix.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Utils/Matrix.hs @@ -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 diff --git a/src/Bio/Utils/Monomer.hs b/src/Bio/Utils/Monomer.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Utils/Monomer.hs @@ -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 diff --git a/test/HandcraftedSpec.hs b/test/HandcraftedSpec.hs new file mode 100644 --- /dev/null +++ b/test/HandcraftedSpec.hs @@ -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 diff --git a/test/JuliaSpec.hs b/test/JuliaSpec.hs new file mode 100644 --- /dev/null +++ b/test/JuliaSpec.hs @@ -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" diff --git a/test/Spec.hs b/test/Spec.hs new file mode 100644 --- /dev/null +++ b/test/Spec.hs @@ -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