diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -2,6 +2,21 @@ ## [Unreleased] +## [0.1.2.3] - 2019-12-12 +### Fixed +- Fixes for .mae pasrser. +- Fixes for instance of `StructureModels` for `Mae`. + +## [0.1.2.2] - 2019-11-27 +### Added +- Instance of `StructureModels` for `Mae`. +### Fixed +- Fix for .mae pasrser. + +## [0.1.2.1] - 2019-11-25 +### Added +- Parser for `MAE`. + ## [0.1.2.0] - 2019-09-03 ### Added - Parser for `FASTA`. diff --git a/cobot-io.cabal b/cobot-io.cabal --- a/cobot-io.cabal +++ b/cobot-io.cabal @@ -4,10 +4,10 @@ -- -- see: https://github.com/sol/hpack -- --- hash: 5649e9d10f32f2f4349434935511b7ca5b206829754168dd6c74817d422c0e91 +-- hash: 77c33361ec08da24af8711ddff5d4e8185bfe4b1fd097bb10c985f15857655f4 name: cobot-io -version: 0.1.2.1 +version: 0.1.2.3 synopsis: Biological data file formats and IO description: Please see the README on GitHub at category: Bio @@ -67,7 +67,7 @@ Paths_cobot_io hs-source-dirs: src - default-extensions: DeriveGeneric DeriveFunctor DeriveFoldable DeriveAnyClass FlexibleInstances InstanceSigs MultiParamTypeClasses RecordWildCards ScopedTypeVariables OverloadedStrings TypeFamilies DataKinds ConstraintKinds TypeOperators TemplateHaskell FlexibleContexts + default-extensions: DeriveGeneric DeriveFunctor DeriveFoldable DeriveAnyClass FlexibleInstances InstanceSigs MultiParamTypeClasses RecordWildCards ScopedTypeVariables OverloadedStrings TypeApplications TypeFamilies DataKinds ConstraintKinds TypeOperators TemplateHaskell FlexibleContexts build-depends: array >=0.5 && <0.6 , attoparsec >=0.10 && <0.14 @@ -98,6 +98,7 @@ GBParserSpec GBWriterSpec MAEParserSpec + MAESpec MMTFSpec SequenceSpec UniprotSpec diff --git a/src/Bio/MAE.hs b/src/Bio/MAE.hs --- a/src/Bio/MAE.hs +++ b/src/Bio/MAE.hs @@ -1,17 +1,38 @@ +{-# LANGUAGE AllowAmbiguousTypes #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} + module Bio.MAE - ( module T + ( Block (..), FromMaeValue (..) + , Mae (..), MaeValue (..) + , Table (..) , fromFile , fromText , maeP ) where import Bio.MAE.Parser -import Bio.MAE.Type as T +import Bio.MAE.Type (Block (..), FromMaeValue (..), + Mae (..), MaeValue (..), Table (..)) +import Bio.Structure (Atom (..), Bond (..), Chain (..), + GlobalID (..), LocalID (..), + Model (..), Residue (..), + SecondaryStructure (..), + StructureModels (..)) import Control.Monad.IO.Class (MonadIO, liftIO) import Data.Attoparsec.Text (parseOnly) -import Data.Bifunctor (first) -import Data.Text (Text, pack) +import Data.Bifunctor (bimap, first) +import Data.Function (on) +import qualified Data.List as L (find, groupBy, sortOn) +import Data.Map.Strict (Map) +import qualified Data.Map.Strict as M (fromList, lookup, member, (!)) +import Data.Maybe (catMaybes, fromJust) +import Data.Text (Text) +import qualified Data.Text as T (head, init, last, null, pack, + strip, tail) import qualified Data.Text.IO as TIO (readFile) +import Data.Vector (Vector) +import qualified Data.Vector as V (fromList) +import Linear.V3 (V3 (..)) -- | Reads 'Mae' from givem file. -- @@ -21,4 +42,131 @@ -- | Reads 'Mae' from 'Text'. -- fromText :: Text -> Either Text Mae -fromText = first pack . parseOnly maeP +fromText = first T.pack . parseOnly maeP + +instance StructureModels Mae where + modelsOf Mae{..} = V.fromList $ fmap blockToModel blocks + where + unsafeGetFromContentsMap :: FromMaeValue a => Map Text [MaeValue] -> Text -> Int -> a + unsafeGetFromContentsMap m name i = unsafeFromMaeValue $ (m M.! name) !! i + + getFromContentsMap :: FromMaeValue a => Map Text [MaeValue] -> Text -> Int -> Maybe a + getFromContentsMap m name i = fromMaeValue $ (m M.! name) !! i + + blockToModel :: Block -> Model + blockToModel Block{..} = Model (atomsTableToChains atomsTable) bonds + where + atomsTable = findTable "m_atom" + numberOfAtoms = length $ atomsTable M.! "r_m_x_coord" + + bondsTable = findTable "m_bond" + (bondGraph, bonds) = bondsTableToGlobalBonds bondsTable + + findTable :: Text -> Map Text [MaeValue] + findTable name = contents $ fromJust $ L.find ((== name) . tableName) tables + + stripQuotes :: Text -> Text + stripQuotes t | not (T.null t) && T.head t == T.last t, T.last t == '\"' = T.strip $ T.init $ T.tail t + | otherwise = T.strip t + + toGroupsOn :: (Eq b, Ord b) => (a -> b) -> [a] -> [[a]] + toGroupsOn f = L.groupBy ((==) `on` f) . L.sortOn f + + bondsTableToGlobalBonds :: Map Text [MaeValue] -> (Map Int [(Int, Int)], Vector (Bond GlobalID)) + bondsTableToGlobalBonds m = bimap toMap V.fromList bonds' + where + numberOfBonds = length $ m M.! "i_m_from" + bonds' = unzip $ fmap indexToBond [0 .. numberOfBonds - 1] + + toMap :: [(Int, (Int, Int))] -> Map Int [(Int, Int)] + toMap = M.fromList . fmap (\l@((k, _) : _) -> (k, fmap snd l)) . toGroupsOn fst + + indexToBond :: Int -> ((Int, (Int, Int)), Bond GlobalID) + indexToBond i = ((x, (y, o)), Bond (GlobalID x) (GlobalID y) o) + where + x = getFromContentsI "i_m_from" - 1 + y = getFromContentsI "i_m_to" - 1 + o = getFromContentsI "i_m_order" + + getFromContentsI :: FromMaeValue a => Text -> a + getFromContentsI = flip (unsafeGetFromContentsMap m) i + + atomsTableToChains :: Map Text [MaeValue] -> Vector Chain + atomsTableToChains m = V.fromList $ fmap groupToChain groupedByChains + where + groupedByChains = toGroupsOn (unsafeGetFromContents @Text "s_m_chain_name") [0 .. numberOfAtoms - 1] + + unsafeGetFromContents :: FromMaeValue a => Text -> Int -> a + unsafeGetFromContents = unsafeGetFromContentsMap m + + getFromContents :: FromMaeValue a => Text -> Int -> Maybe a + getFromContents = getFromContentsMap m + + groupToChain :: [Int] -> Chain + groupToChain [] = error "Group that is result of List.groupBy can't be empty." + groupToChain group@(h : _) = Chain name residues + where + name = stripQuotes $ unsafeGetFromContents "s_m_chain_name" h + + groupedByResidues = toGroupsOn by group + residues = V.fromList $ fmap groupToResidue groupedByResidues + + by :: Int -> (Int, Text) + by i = (unsafeGetFromContents "i_m_residue_number" i, unsafeGetFromContents "s_m_pdb_residue_name" i) + + groupToResidue :: [Int] -> Residue + groupToResidue [] = error "Group that is result of List.groupBy can't be empty." + groupToResidue group@(h : _) = Residue name atoms (V.fromList localBonds) secondary chemCompType + where + name = stripQuotes $ unsafeGetFromContents "s_m_pdb_residue_name" h + atoms = V.fromList $ fmap indexToAtom group + + localInds = [0 .. length group - 1] + globalToLocal = M.fromList $ zip group localInds + bondsParts = fmap (`M.lookup` bondGraph) group + localBonds = concat $ catMaybes $ zipWith (\l x -> fmap (concatMap (toLocalBond x)) l) bondsParts localInds + + toLocalBond :: Int -> (Int, Int) -> [Bond LocalID] + toLocalBond x (y, o) | y `elem` group = pure $ Bond (LocalID x) + (LocalID $ globalToLocal M.! y) + o + | otherwise = [] + + secondary = Undefined + chemCompType = mempty + + indexToAtom :: Int -> Atom + indexToAtom i = Atom (GlobalID i) + (stripQuotes $ getFromContentsI "s_m_pdb_atom_name") + (elIndToElement M.! getFromContentsI "i_m_atomic_number") + coords + (getFromContentsIWithDef 0 "i_m_formal_charge") + (getFromContentsIWithDef 0 "r_m_pdb_tfactor") + (getFromContentsIWithDef 0 "r_m_pdb_occupancy") + where + getFromContentsI :: FromMaeValue a => Text -> a + getFromContentsI = flip unsafeGetFromContents i + + getFromContentsIWithDef :: FromMaeValue a => a -> Text -> a + getFromContentsIWithDef def n | n `M.member` m = maybe def id $ getFromContents n i + | otherwise = def + + coords :: V3 Float + coords = V3 (getFromContentsI "r_m_x_coord") + (getFromContentsI "r_m_y_coord") + (getFromContentsI "r_m_z_coord") + +elIndToElement :: Map Int Text +elIndToElement = M.fromList $ zip [1 .. 118] [ "H", "He", "Li", "Be", "B", "C", "N", "O", "F", "Ne" + , "Na", "Mg", "Al", "Si", "P", "S", "Cl", "Ar", "K", "Ca" + , "Sc", "Ti", "V", "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn" + , "Ga", "Ge", "As", "Se", "Br", "Kr", "Rb", "Sr", "Y", "Zr" + , "Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", "Cd", "In", "Sn" + , "Sb", "Te", "I", "Xe", "Cs", "Ba", "La", "Ce", "Pr", "Nd" + , "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb" + , "Lu", "Hf", "Ta", "W", "Re", "Os", "Ir", "Pt", "Au", "Hg" + , "Tl", "Pb", "Bi", "Po", "At", "Rn", "Fr", "Ra", "Ac", "Th" + , "Pa", "U", "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm" + , "Md", "No", "Lr", "Rf", "Db", "Sg", "Bh", "Hs", "Mt", "Ds" + , "Rg", "Cn", "Nh", "Fl", "Mc", "Lv", "Ts", "Og" + ] diff --git a/src/Bio/MAE/Parser.hs b/src/Bio/MAE/Parser.hs --- a/src/Bio/MAE/Parser.hs +++ b/src/Bio/MAE/Parser.hs @@ -11,15 +11,15 @@ Table (..)) import Control.Applicative ((<|>)) import Control.Monad (replicateM, when, zipWithM) -import Data.Attoparsec.Text (Parser, char, decimal, endOfInput, - endOfLine, many', many1', string, - takeWhile, takeWhile1) +import Data.Attoparsec.Text (Parser, anyChar, char, decimal, + endOfInput, endOfLine, many', many1', + peekChar, string, takeWhile, takeWhile1) import Data.Char (isSpace) import Data.List (transpose) import Data.Map.Strict (Map) import qualified Data.Map.Strict as M (fromList) import Data.Text (Text) -import qualified Data.Text as T (pack, uncons) +import qualified Data.Text as T (cons, pack, uncons) import qualified Data.Text.Read as TR (decimal, rational, signed) import Prelude hiding (takeWhile) @@ -29,9 +29,7 @@ <* endOfInput versionP :: Parser Text -versionP = inBrackets $ lineP - *> delimiterP - *> lineP +versionP = many' tillEndOfLine *> inBrackets (lineP *> delimiterP *> lineP) blockP :: Parser Block blockP = uncurry <$> (Block <$> anyStringP <* many' oneSpaceP) @@ -137,7 +135,7 @@ anyStringP = takeWhile1 (not . isSpace) valueTP :: Parser Text -valueTP = ((<>) <$> string quoteT <*> ((<>) <$> takeWhile (/= quote) <*> string quoteT)) +valueTP = ((<>) <$> string quoteT <*> ((<>) <$> notQuote <*> string quoteT)) <|> anyStringP where quote :: Char @@ -145,6 +143,20 @@ quoteT :: Text quoteT = T.pack $ pure quote + + notQuote :: Parser Text + notQuote = do + curCharPeek <- peekChar + case curCharPeek of + Just '\\' -> do + curChar <- anyChar + nextCharPeek <- peekChar + case nextCharPeek of + Just '\"' -> anyChar >>= \x -> fmap (T.cons curChar . T.cons x) notQuote + _ -> notQuote >>= pure . T.cons curChar + Just '\"' -> pure mempty + Just _ -> anyChar >>= \x -> fmap (T.cons x) notQuote + Nothing -> pure mempty commentaryP :: Parser () commentaryP = () <$ many' (many' oneSpaceP *> char '#' *> takeWhile (`notElem` ['\n', '\r']) *> endOfLine) diff --git a/test/MAEParserSpec.hs b/test/MAEParserSpec.hs --- a/test/MAEParserSpec.hs +++ b/test/MAEParserSpec.hs @@ -43,6 +43,9 @@ it "h2o.mae" $ parseFileSpec "test/MAE/h2o.mae" it "docking_1.mae" $ parseFileSpec "test/MAE/docking_1.mae" it "docking_2.mae" $ parseFileSpec "test/MAE/docking_2.mae" + it "docking_1_2.mae" $ parseFileSpec "test/MAE/docking_1_2.mae" + it "schrod_test.mae" $ parseFileSpec "test/MAE/schrod_test.mae" + it "GLY_3.mae" $ parseFileSpec "test/MAE/GLY_3.mae" parseFileSpec :: FilePath -> Expectation parseFileSpec path = fromFile path >> pure () @@ -89,7 +92,7 @@ tableWithQuotedMaeValuesMap :: [(Text, [MaeValue])] tableWithQuotedMaeValuesMap = [ ("s_quote0", fmap StringMaeValue ["\" ssss s \"", "\" more quotes \"", "\" this is long and dull comment \""]) , ("r_neg", fmap RealMaeValue [-1.0, -2.28, -3.22]) - , ("s_quote", fmap StringMaeValue ["\" \"", "aaa", "\" this is long and dull comment \""]) + , ("s_quote", fmap StringMaeValue ["\" \"", "aaa", "\" this is long \\\" and dull comment \""]) , ("r_neg", fmap RealMaeValue [-1.0, -2.28, -3.22]) ] @@ -97,7 +100,7 @@ tableWithQuotedMaeValues = Table "table_with_quoted_values" $ M.fromList tableWithQuotedMaeValuesMap tableWithQuotedMaeValuesT :: Text -tableWithQuotedMaeValuesT = "table_with_quoted_values[3]{ \n s_quote0\n r_neg \n s_quote\n r_neg \n ::: \n 1 \" ssss s \" -1.0 \" \" -1.0 \n 2 \" more quotes \" -2.28 aaa -2.28 \n 3 \" this is long and dull comment \" -3.22 \" this is long and dull comment \" -3.22\n ::: \n} \n" +tableWithQuotedMaeValuesT = "table_with_quoted_values[3]{ \n s_quote0\n r_neg \n s_quote\n r_neg \n ::: \n 1 \" ssss s \" -1.0 \" \" -1.0 \n 2 \" more quotes \" -2.28 aaa -2.28 \n 3 \" this is long and dull comment \" -3.22 \" this is long \\\" and dull comment \"-3.22\n ::: \n} \n" tableWithCommentsMap :: [(Text, [MaeValue])] tableWithCommentsMap = [ ("i_val", replicate 7 Absent) diff --git a/test/MAESpec.hs b/test/MAESpec.hs new file mode 100644 --- /dev/null +++ b/test/MAESpec.hs @@ -0,0 +1,98 @@ +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE RecordWildCards #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} + +module MAESpec where + +import Bio.MAE (fromFile) +import Bio.Structure (Atom (..), Bond (..), Chain (..), GlobalID (..), + LocalID (..), Model (..), Residue (..), + StructureModels (..)) +import Data.Set (Set) +import qualified Data.Set as S (fromList) +import Data.Vector (Vector) +import qualified Data.Vector as V (fromList, toList, (!)) +import Linear.V3 (V3 (..)) +import Test.Hspec + +maeSpec :: Spec +maeSpec = describe "Mae spec." $ do + Model{..} <- runIO $ V.toList . modelsOf <$> fromFile "test/MAE/small.mae" >>= \[x] -> pure x + + it "two chains" $ length modelChains `shouldBe` 2 + it "residues" $ do + fmap resName (chainResidues $ modelChains V.! 0) `shouldBe` V.fromList ["ACE", "ASP", "ILE", "LYS"] + fmap resName (chainResidues $ modelChains V.! 1) `shouldBe` V.fromList ["GLU", "LEU", "VAL", "ARG", "PRO", "GLY", "ALA", "LEU", "VAL"] + + it "atoms count" $ do + sum (fmap (length . resAtoms) $ chainResidues $ modelChains V.! 0) `shouldBe` 62 + sum (fmap (length . resAtoms) $ chainResidues $ modelChains V.! 1) `shouldBe` 140 + + let allBonds = [ (2, 15, 1), (31, 44, 1), (15, 60, 1), (52, 27, 1), (57, 30, 1), (56, 8, 1), (38, 9, 1), (16, 36, 1), (31, 35, 1), (42, 22, 1), (38, 10, 1), (46, 40, 1), (36, 8, 1) + , (46, 7, 1), (35, 13, 1), (60, 31, 1), (40, 32, 1), (46, 45, 1), (29, 1, 1), (9, 53, 1), (22, 33, 1), (41, 18, 1), (45, 31, 1), (13, 53, 1), (40, 19, 1), (55, 20, 1) + , (16, 22, 1), (51, 36, 1), (46, 56, 1), (8, 17, 1), (56, 40, 1), (50, 6, 1), (57, 16, 1), (57, 34, 1), (48, 47, 1), (28, 48, 1), (57, 30, 1), (39, 40, 1), (11, 6, 1) + , (35, 55, 1), (47, 30, 1), (4, 35, 1), (60, 35, 1), (39, 56, 1), (44, 24, 1), (29, 55, 1), (29, 41, 1), (6, 55, 1), (52, 51, 1), (32, 21, 1), (55, 55, 1), (43, 55, 1) + , (30, 44, 1), (54, 47, 1), (13, 50, 1), (14, 56, 1), (44, 54, 1), (58, 21, 1), (27, 58, 1), (84, 122, 2), (122, 122, 2), (109, 79, 2), (103, 121, 2), (102, 139, 2) + , (93, 93, 2), (68, 103, 2), (110, 77, 2), (109, 63, 2), (86, 83, 2), (123, 137, 2), (110, 131, 2), (122, 85, 2), (75, 110, 2), (131, 138, 2), (88, 134, 2), (117, 81, 2) + , (107, 74, 2), (67, 61, 2), (97, 134, 2), (94, 131, 2), (65, 95, 2), (124, 100, 2), (120, 106, 2), (71, 111, 2), (95, 129, 2), (104, 116, 2), (61, 95, 2), (88, 132, 2) + , (80, 97, 2), (68, 121, 2), (138, 87, 2), (84, 134, 2), (86, 139, 2), (84, 76, 2), (83, 85, 2), (76, 61, 2), (116, 98, 2), (64, 82, 2), (106, 93, 2), (96, 102, 2), (98, 122, 2) + , (74, 82, 2), (123, 130, 2), (114, 127, 2), (102, 122, 2), (90, 126, 2), (118, 92, 2), (88, 71, 2), (120, 132, 2), (135, 71, 2), (125, 136, 2), (71, 67, 2), (128, 110, 2) + , (67, 95, 2), (83, 72, 2), (60, 70, 2), (99, 135, 2), (87, 87, 2), (63, 124, 2), (64, 110, 2), (127, 81, 2), (90, 106, 2), (93, 92, 2), (62, 91, 2), (67, 119, 2), (95, 70, 2) + , (134, 115, 2), (113, 115, 2), (93, 112, 2), (102, 106, 2), (136, 105, 2), (82, 126, 2), (109, 103, 2), (66, 119, 2), (84, 60, 2), (79, 67, 2), (63, 96, 2), (74, 134, 2), (112, 91, 2) + , (105, 68, 2), (125, 80, 2), (136, 75, 2), (62, 117, 2), (66, 124, 2), (105, 120, 2), (141, 63, 2), (63, 90, 2), (139, 60, 2), (77, 73, 2), (81, 118, 2), (81, 72, 2), (123, 126, 2) + , (128, 136, 2), (83, 91, 2), (86, 124, 2), (140, 108, 2), (126, 98, 2), (67, 92, 2), (131, 69, 2), (69, 94, 2), (137, 66, 2), (99, 87, 2), (134, 134, 2), (134, 72, 2), (91, 138, 2) + , (89, 97, 2), (109, 128, 2), (123, 100, 2), (61, 64, 2), (68, 114, 2), (119, 115, 2), (122, 106, 2), (131, 84, 2), (93, 94, 2), (124, 118, 2), (82, 119, 2), (107, 95, 2), (112, 64, 2) + , (71, 71, 2), (89, 62, 2), (89, 61, 2), (120, 80, 2), (85, 61, 2), (97, 95, 2), (94, 98, 2), (98, 139, 2), (115, 133, 2), (115, 64, 2), (136, 66, 2), (131, 62, 2), (83, 102, 2) + , (121, 60, 2), (117, 110, 2), (134, 138, 2), (119, 89, 2), (137, 128, 2), (112, 118, 2), (105, 116, 2), (95, 69, 2) + ] + + it "global bonds" $ do + modelBonds `shouldBe` toBond GlobalID <$> V.fromList allBonds + + let residue1 = (chainResidues $ modelChains V.! 0) V.! 3 + let residue2 = (chainResidues $ modelChains V.! 1) V.! 0 + + it "atoms in residues" $ do + let atoms1 = resAtoms residue1 + let atoms2 = resAtoms residue2 + + fmap atomId atoms1 `shouldBe` GlobalID <$> V.fromList ([37 .. 58] <> [199 .. 201]) + fmap atomId atoms2 `shouldBe` GlobalID <$> V.fromList [59 .. 73] + + fmap atomName atoms1 `shouldBe` V.fromList [ "N", "CA", "C", "O", "CB", "CG", "CD", "CE", "NZ", "H" + , "HA", "HB3", "HB2", "HG3", "HG2", "HD3", "HD2", "HE3" + , "HE2", "HZ1", "HZ2", "HZ3", "2H", "CG2", "HA" + ] + fmap atomName atoms2 `shouldBe` V.fromList [ "N", "CA", "C", "O", "CB", "CG", "CD", "OE1", "OE2", "H" + , "HA", "HB3", "HB2", "HG3", "HG2" + ] + + fmap atomCoords atoms1 `shouldBe` V.fromList [ V3 (-3.444000) 7.750000 48.589000, V3 (-2.336000) 8.349000 47.850000, V3 (-1.127000) 7.408000 47.751000, V3 (-0.731000) 6.818000 48.755000, V3 (-1.974000) 9.706000 48.503000 + , V3 (-0.665000) 10.349000 47.988000, V3 (-0.475000) 11.818000 48.390000, V3 (-1.315000) 12.784000 47.541000, V3 (-1.050000) 14.186000 47.903000, V3 (-3.565000) 8.072000 49.540000 + , V3 (-2.692000) 8.558000 46.839000, V3 (-1.882000) 9.581000 49.583000, V3 (-2.810000) 10.388000 48.353000, V3 (-0.610000) 10.265000 46.901000, V3 0.184000 9.784000 48.375000 + , V3 0.582000 12.073000 48.302000, V3 (-0.729000) 11.939000 49.444000, V3 (-2.380000) 12.585000 47.665000, V3 (-1.084000) 12.653000 46.483000, V3 (-1.609000) 14.795000 47.322000 + , V3 (-1.289000) 14.335000 48.873000, V3 (-0.072000) 14.393000 47.762000, V3 (-9.600000) 7.518000 44.746000, V3 (-6.594000) 4.327000 47.986000, V3 (-5.134000) 6.683000 50.072000 + ] + fmap atomCoords atoms2 `shouldBe` V.fromList [ V3 (-8.382000) 11.633000 16.946000, V3 (-9.715000) 12.157000 17.191000, V3 (-9.590000) 13.665000 17.450000, V3 (-9.030000) 14.054000 18.475000, V3 (-10.323000) 11.388000 18.396000 + , V3 (-11.833000) 11.576000 18.604000, V3 (-12.658000) 10.923000 17.495000, V3 (-12.525000) 9.689000 17.338000, V3 (-13.405000) 11.669000 16.827000, V3 (-7.615000) 12.187000 17.300000 + , V3 (-10.322000) 11.996000 16.297000, V3 (-9.798000) 11.676000 19.305000, V3 (-10.128000) 10.318000 18.313000, V3 (-12.084000) 12.633000 18.684000, V3 (-12.123000) 11.111000 19.544000 + ] + + it "bonds in residues" $ do + toSet (resBonds residue1) `shouldBe` toSet (V.fromList $ toBond LocalID . (\(x, y, o) -> (x - 37, y - 37, o)) <$> filter (\(x, y, _) -> x - 1 `elem` [37 .. 58] && y - 1 `elem` [37 .. 58]) allBonds) + toSet (resBonds residue2) `shouldBe` toSet (V.fromList $ toBond LocalID . (\(x, y, o) -> (x - 59, y - 59, o)) <$> filter (\(x, y, _) -> x - 1 `elem` [59 .. 73] && y - 1 `elem` [59 .. 73]) allBonds) + where + toBond :: (Int -> b) -> (Int, Int, Int) -> Bond b + toBond f (x, y, o) = Bond (f $ x - 1) (f $ y - 1) o + + toSet :: Ord a => Vector a -> Set a + toSet = S.fromList . V.toList + +instance Ord (Bond LocalID) where + (Bond (LocalID x) (LocalID y) _) <= (Bond (LocalID x') (LocalID y') _) | x == x' = y <= y' + | otherwise = x <= x' + +instance Ord (Bond GlobalID) where + (Bond (GlobalID x) (GlobalID y) _) <= (Bond (GlobalID x') (GlobalID y') _) | x == x' = y <= y' + | otherwise = x <= x' diff --git a/test/Spec.hs b/test/Spec.hs --- a/test/Spec.hs +++ b/test/Spec.hs @@ -5,6 +5,7 @@ import GBParserSpec import GBWriterSpec import MAEParserSpec +import MAESpec import MMTFSpec import SequenceSpec import System.IO @@ -38,3 +39,4 @@ fastaWriterSpec -- Mae maeParserSpec + maeSpec