diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -2,6 +2,13 @@ ## [Unreleased] +## [0.1.3.3] - 2020-05-07 +### Added +- Type and decoder for `ab1` with raw channel data and peak locations. +### Fixed +- Use supplied thresholds instead of default ones in the implementation of `cleanWith`; +- Do not run `IO` in `hspec` `Spec` monad, use `beforeAll`. + ## [0.1.3.2] - 2020-04-15 ### Changed - Resolver version up. diff --git a/cobot-io.cabal b/cobot-io.cabal --- a/cobot-io.cabal +++ b/cobot-io.cabal @@ -1,136 +1,140 @@ cabal-version: 1.12 -name: cobot-io -version: 0.1.3.2 -license: BSD3 -license-file: LICENSE -copyright: 2020, Biocad -maintainer: pavel@yakovlev.me -author: Pavel Yakovlev, Bogdan Neterebskii, Alexander Sadovnikov -homepage: https://github.com/biocad/cobot-io#readme -bug-reports: https://github.com/biocad/cobot-io/issues -synopsis: Biological data file formats and IO -description: - Please see the README on GitHub at -category: Bio -build-type: Simple + +-- This file has been generated from package.yaml by hpack version 0.33.0. +-- +-- see: https://github.com/sol/hpack +-- +-- hash: 476c2b78718477cf61d3cc8b8f752cbb5b4533d9a349a7aa05dc96fa352f2d7e + +name: cobot-io +version: 0.1.3.3 +synopsis: Biological data file formats and IO +description: Please see the README on GitHub at +category: Bio +homepage: https://github.com/biocad/cobot-io#readme +bug-reports: https://github.com/biocad/cobot-io/issues +author: Pavel Yakovlev, Bogdan Neterebskii, Alexander Sadovnikov +maintainer: pavel@yakovlev.me +copyright: 2020, Biocad +license: BSD3 +license-file: LICENSE +build-type: Simple extra-source-files: README.md ChangeLog.md source-repository head - type: git - location: https://github.com/biocad/cobot-io + type: git + location: https://github.com/biocad/cobot-io library - exposed-modules: - Bio.ABI - Bio.ABI.Clean - Bio.ABI.Decode - Bio.FASTA - Bio.FASTA.Parser - Bio.FASTA.Type - Bio.FASTA.Writer - Bio.GB - Bio.GB.Parser - Bio.GB.Type - Bio.GB.Writer - Bio.MAE - Bio.MAE.Parser - Bio.MAE.Type - Bio.MMTF - Bio.MMTF.Decode - Bio.MMTF.Decode.Codec - Bio.MMTF.Decode.MessagePack - Bio.MMTF.MessagePack - Bio.MMTF.Type - Bio.PDB - Bio.PDB.BondRestoring - Bio.PDB.Functions - Bio.PDB.Parser - Bio.PDB.Reader - Bio.PDB.Type - Bio.Sequence - Bio.Sequence.Basecalled - Bio.Sequence.Basecalled.Type - Bio.Sequence.Class - Bio.Sequence.Functions.Marking - Bio.Sequence.Functions.Sequence - Bio.Sequence.Functions.Weight - Bio.Sequence.Utilities - Bio.Structure - Bio.Structure.Functions - Bio.Uniprot - Bio.Uniprot.Parser - Bio.Uniprot.Type - hs-source-dirs: src - other-modules: - Paths_cobot_io - default-language: Haskell2010 - default-extensions: DeriveGeneric DeriveFunctor DeriveFoldable - DeriveAnyClass FlexibleInstances InstanceSigs MultiParamTypeClasses - RecordWildCards ScopedTypeVariables OverloadedStrings - TypeApplications TypeFamilies DataKinds ConstraintKinds - TypeOperators TemplateHaskell FlexibleContexts - build-depends: - array ==0.5.*, - attoparsec >=0.10 && <0.14, - base >=4.7 && <5, - binary >=0.8.3.0 && <1.0, - bytestring >=0.10.8.1 && <0.11, - containers >=0.5.7.1 && <0.7, - data-msgpack >=0.0.9 && <0.1, - deepseq ==1.4.*, - http-conduit ==2.3.*, - hyraxAbif >=0.2.3.15 && <0.2.4.0, - lens >=4.16 && <5.0, - linear ==1.20.*, - mtl >=2.2.1 && <2.3.0, - split >=0.2.3.4 && <0.3, - text >=1.2.2.1 && <1.3, - vector >=0.12.1.2 && <0.13 + exposed-modules: + Bio.ABI + Bio.ABI.Clean + Bio.ABI.Decode + Bio.FASTA + Bio.FASTA.Parser + Bio.FASTA.Type + Bio.FASTA.Writer + Bio.GB + Bio.GB.Parser + Bio.GB.Type + Bio.GB.Writer + Bio.MAE + Bio.MAE.Parser + Bio.MAE.Type + Bio.MMTF + Bio.MMTF.Decode + Bio.MMTF.Decode.Codec + Bio.MMTF.Decode.MessagePack + Bio.MMTF.MessagePack + Bio.MMTF.Type + Bio.PDB + Bio.PDB.BondRestoring + Bio.PDB.Functions + Bio.PDB.Parser + Bio.PDB.Reader + Bio.PDB.Type + Bio.Sequence + Bio.Sequence.Basecalled + Bio.Sequence.Basecalled.Type + Bio.Sequence.Class + Bio.Sequence.Functions.Marking + Bio.Sequence.Functions.Sequence + Bio.Sequence.Functions.Weight + Bio.Sequence.Utilities + Bio.Structure + Bio.Structure.Functions + Bio.Uniprot + Bio.Uniprot.Parser + Bio.Uniprot.Type + other-modules: + Paths_cobot_io + hs-source-dirs: + src + 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 + , base >=4.7 && <5 + , binary >=0.8.3.0 && <1.0 + , bytestring >=0.10.8.1 && <0.11 + , containers >=0.5.7.1 && <0.7 + , data-msgpack >=0.0.9 && <0.1 + , deepseq >=1.4 && <1.5 + , http-conduit >=2.3 && <2.4 + , hyraxAbif >=0.2.3.15 && <0.2.4.0 + , lens >=4.16 && <5.0 + , linear >=1.20 && <1.21 + , mtl >=2.2.1 && <2.3.0 + , split + , text >=1.2.2.1 && <1.3 + , vector + default-language: Haskell2010 test-suite cobot-io-test - type: exitcode-stdio-1.0 - main-is: Spec.hs - hs-source-dirs: test - other-modules: - ABISpec - FastaParserSpec - FASTASpec - FastaWriterSpec - GBParserSpec - GBWriterSpec - MAEParserSpec - MAESpec - MMTFSpec - PDBParserSpec - PDBSpec - SequenceSpec - StructureSpec - UniprotSpec - Paths_cobot_io - default-language: Haskell2010 - default-extensions: OverloadedStrings TypeFamilies - ghc-options: -threaded -rtsopts -with-rtsopts=-N - build-depends: - QuickCheck >=2.9.2 && <2.14, - array ==0.5.*, - attoparsec >=0.10 && <0.14, - base >=4.7 && <5, - binary >=0.8.3.0 && <1.0, - bytestring >=0.10.8.1 && <0.11, - cobot-io -any, - containers >=0.5.7.1 && <0.7, - data-msgpack >=0.0.9 && <0.1, - deepseq ==1.4.*, - directory >=1.3.6.0 && <1.4, - hspec >=2.4.1 && <2.8, - http-conduit ==2.3.*, - hyraxAbif >=0.2.3.15 && <0.2.4.0, - lens >=4.16 && <5.0, - linear ==1.20.*, - mtl >=2.2.1 && <2.3.0, - neat-interpolation ==0.3.*, - split >=0.2.3.4 && <0.3, - text >=1.2.2.1 && <1.3, - vector >=0.12.1.2 && <0.13 + type: exitcode-stdio-1.0 + main-is: Spec.hs + other-modules: + ABISpec + FastaParserSpec + FASTASpec + FastaWriterSpec + GBParserSpec + GBWriterSpec + MAEParserSpec + MAESpec + MMTFSpec + PDBParserSpec + PDBSpec + SequenceSpec + StructureSpec + UniprotSpec + Paths_cobot_io + hs-source-dirs: + test + default-extensions: OverloadedStrings TypeFamilies + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: + QuickCheck >=2.9.2 && <2.14 + , array >=0.5 && <0.6 + , attoparsec >=0.10 && <0.14 + , base >=4.7 && <5 + , binary >=0.8.3.0 && <1.0 + , bytestring >=0.10.8.1 && <0.11 + , cobot-io + , containers >=0.5.7.1 && <0.7 + , data-msgpack >=0.0.9 && <0.1 + , deepseq >=1.4 && <1.5 + , directory + , hspec >=2.4.1 && <2.8 + , http-conduit >=2.3 && <2.4 + , hyraxAbif >=0.2.3.15 && <0.2.4.0 + , lens >=4.16 && <5.0 + , linear >=1.20 && <1.21 + , mtl >=2.2.1 && <2.3.0 + , neat-interpolation >=0.3 + , split + , text >=1.2.2.1 && <1.3 + , vector + default-language: Haskell2010 diff --git a/src/Bio/ABI/Clean.hs b/src/Bio/ABI/Clean.hs --- a/src/Bio/ABI/Clean.hs +++ b/src/Bio/ABI/Clean.hs @@ -1,6 +1,8 @@ module Bio.ABI.Clean ( Cleanable (..) + , Thresholds (..) + , defaultThresholds ) where import Bio.Sequence (mean, meanInRange) @@ -18,22 +20,25 @@ cleanWith :: Thresholds -> a -> Maybe a -- | Thresholds to clean the data. +-- -- ABI file contains sequence with quality. -- By design of sanger sequencing method start and end of the sequence have bad quality. -- Moreover, internal part of the sequence can also has bad quality. -- To clean the data we make 2 steps. -- -- Step 1. Clean edges: +-- -- * take frame with @frameSize@ and go through the sequence; -- * on each step evaluate mean value; --- * if mean value less than @edgeThreshold@, go further; --- * if mean value more than @edgeThreshold@, stop and cut the sequence from END of this frame; +-- * if mean value less than 'edgeThreshold', go further; +-- * if mean value more than 'edgeThreshold' stop and cut the sequence from END of this frame; -- * repeat this algorithm for the right edge. -- -- Step 2. Evaluate quality: +-- -- * for cropped sequence evaluate mean value; --- * if mean value less then @innerThreshold@, sequence is bad; --- * if mean value more then @innerThreshold@, sequence is acceptable. +-- * if mean value less then 'innerThreshold', sequence is bad; +-- * if mean value more then 'innerThreshold', sequence is acceptable. -- -- Logic of this algorithm and 'defaultThresholds' were obtained by taking experiments with read ABI files. -- @@ -54,10 +59,10 @@ then fromBoth else Nothing where - fromLeft = cutEdge defaultThresholds input + fromLeft = cutEdge thr input fromBoth = fmap S.reverse . join - $ cutEdge defaultThresholds + $ cutEdge thr . S.reverse <$> fromLeft diff --git a/src/Bio/ABI/Decode.hs b/src/Bio/ABI/Decode.hs --- a/src/Bio/ABI/Decode.hs +++ b/src/Bio/ABI/Decode.hs @@ -1,19 +1,55 @@ {-# OPTIONS_GHC -fno-warn-orphans #-} -module Bio.ABI.Decode () where +module Bio.ABI.Decode + ( decodeRawSequence + , decodeRawSequence' + ) where -import Bio.Sequence (SequenceDecodable (..), - weightedSequence) -import Bio.Sequence.Basecalled (BasecalledSequence) +import Control.Applicative (many) +import Data.Bifunctor (bimap) +import Data.Binary.Get (getInt16be, runGetOrFail) import Data.ByteString as BS (ByteString) import Data.ByteString.Lazy as BSL (ByteString, fromStrict) import Data.ByteString.Lazy.Char8 as BSL8 (unpack) import Data.Char (ord) +import Data.Functor ((<&>)) +import Data.Int (Int16) import Data.List (find) -import Data.Text (Text) +import Data.Text (Text, pack) +import Data.Vector (Vector, fromList) import Hyrax.Abif (Abif (..), Directory (..)) import Hyrax.Abif.Read (getAbif) +import Bio.Sequence (SequenceDecodable (..), + weightedSequence) +import Bio.Sequence.Basecalled (BasecalledSequence, + BasecalledSequenceWithRawData (..)) + +-- | Decode ABIF file with additional raw data that may be required for later processing. +decodeRawSequence :: BSL.ByteString -> Either Text BasecalledSequenceWithRawData +decodeRawSequence bs = do + abif <- getAbif bs + sequence' <- extractSequence abif + quality' <- extractQuality abif + bsSequence <- weightedSequence sequence' quality' + + bsRawG <- findDataByDirectory "DATA" 9 abif >>= decodeShortArray + bsRawA <- findDataByDirectory "DATA" 10 abif >>= decodeShortArray + bsRawT <- findDataByDirectory "DATA" 11 abif >>= decodeShortArray + bsRawC <- findDataByDirectory "DATA" 12 abif >>= decodeShortArray + + bsPeakLocations <- findDataByDirectory "PLOC" 2 abif >>= decodeShortArray <&> fmap fromIntegral + + return BasecalledSequenceWithRawData{..} + +-- | Same as 'decodeRawSequence', for strict @ByteString@. +decodeRawSequence' :: BS.ByteString -> Either Text BasecalledSequenceWithRawData +decodeRawSequence' = decodeRawSequence . BSL.fromStrict + +-- | Discards raw data information. +instance SequenceDecodable BasecalledSequenceWithRawData BasecalledSequence where + sequenceDecode = pure . bsSequence + -- | Converts 'Data.ByteString.Lazy.ByteString' (that should be content of ABI file) -- into 'BasecalledSequence'. -- @@ -39,13 +75,13 @@ -- | Extracts sequence from ABI file. -- extractSequence :: Abif -> Either Text String -extractSequence abif = findDataByDirectory "PBAS" abif >>= checkACGT +extractSequence abif = findDataByDirectory "PBAS" 1 abif <&> BSL8.unpack . dData >>= checkACGT -- | Extracts quality from ABI file. -- Number are encoded with letters, thus we have function @fromIntegral . ord@. -- extractQuality :: Abif -> Either Text [Double] -extractQuality abif = map (fromIntegral . ord) <$> findDataByDirectory "PCON" abif +extractQuality abif = map (fromIntegral . ord) . BSL8.unpack . dData <$> findDataByDirectory "PCON" 1 abif -- | Checks that all chars are from alphabet ACGT -- @@ -58,13 +94,22 @@ -- | Looks into ABI file and extract data by 'Directory' name. -- -findDataByDirectory :: Text -> Abif -> Either Text String -findDataByDirectory dirName abif = - let directoryM = find (\Directory{..} -> dTagName == dirName) . aDirs $ abif - in maybe (Left errorMsg) (Right . getData) directoryM +findDataByDirectory + :: Text -- ^ Directory name + -> Int -- ^ 1-based directory index + -> Abif + -> Either Text Directory +findDataByDirectory dirName dirIndex abif = + let directoryM = find (\Directory{..} -> dTagName == dirName && dTagNum == dirIndex) . aDirs $ abif + in maybe (Left errorMsg) Right directoryM where errorMsg :: Text - errorMsg = "Bio.ABI.Decode: could not find directory " <> dirName + errorMsg = "Bio.ABI.Decode: could not find directory " <> dirName <> " with index " <> pack (show dirIndex) - getData :: Directory -> String - getData = BSL8.unpack . dData +decodeShortArray :: Directory -> Either Text (Vector Int16) +decodeShortArray = + bimap + (\(_, _, msg) -> "Data.ABI.Decode: could not decode short array: " <> pack msg) + (\(_, _, lst) -> fromList lst) + . runGetOrFail (many getInt16be) + . dData diff --git a/src/Bio/Sequence/Basecalled/Type.hs b/src/Bio/Sequence/Basecalled/Type.hs --- a/src/Bio/Sequence/Basecalled/Type.hs +++ b/src/Bio/Sequence/Basecalled/Type.hs @@ -1,13 +1,32 @@ module Bio.Sequence.Basecalled.Type ( BasecalledSequence (..) + , BasecalledSequenceWithRawData (..) ) where import Bio.Sequence (IsSequence (..), WeightedSequence) import Data.Coerce (coerce) +import Data.Int (Int16) +import Data.Vector (Vector) newtype BasecalledSequence = BasecalledSequence (WeightedSequence Double Char) deriving (Eq, Show) +data BasecalledSequenceWithRawData + = BasecalledSequenceWithRawData + { bsSequence :: BasecalledSequence + -- A sequence with quality. + , bsRawG :: Vector Int16 + , bsRawA :: Vector Int16 + , bsRawT :: Vector Int16 + , bsRawC :: Vector Int16 + , bsPeakLocations :: Vector Int + -- ^ Same length as 'bsSequence'. + -- + -- For every base in 'bsSequence' corresponding element in 'bsPeakLocations' is + -- an index in one of the raw vectors. + } + deriving (Eq, Show) + instance IsSequence BasecalledSequence where type Element BasecalledSequence = Char type Marking BasecalledSequence = () @@ -15,4 +34,3 @@ toSequence = coerce fromSequence = coerce - diff --git a/test/ABISpec.hs b/test/ABISpec.hs --- a/test/ABISpec.hs +++ b/test/ABISpec.hs @@ -1,13 +1,18 @@ +{-# LANGUAGE RecordWildCards #-} + module ABISpec where -import Bio.ABI (Cleanable (..)) -import Bio.Sequence (SequenceDecodable (..)) -import qualified Bio.Sequence as S (getWeights, length, toList) -import Bio.Sequence.Basecalled (BasecalledSequence) import Data.ByteString.Lazy as BSL (readFile) import Data.Text (Text) import Test.Hspec +import Bio.ABI (Cleanable (..)) +import Bio.ABI.Decode (decodeRawSequence) +import Bio.Sequence (SequenceDecodable (..)) +import qualified Bio.Sequence as S (getWeights, length, toList) +import Bio.Sequence.Basecalled (BasecalledSequence, + BasecalledSequenceWithRawData (..)) + abiExtractSpec :: Spec abiExtractSpec = describe "ABI decode" $ do @@ -21,6 +26,13 @@ datM <- readData "test/ABI/not_ab1.txt" datM `shouldBe` Left "Error reading root: not enough bytes" + it "decode with raw data" $ do + Right BasecalledSequenceWithRawData{..} <- decodeRawSequence <$> BSL.readFile "test/ABI/test.ab1" + S.length bsSequence `shouldBe` length bsPeakLocations + length bsRawG `shouldSatisfy` (>0) + length bsRawA `shouldSatisfy` (>0) + length bsRawT `shouldSatisfy` (>0) + length bsRawC `shouldSatisfy` (>0) abiCleanSpec :: Spec abiCleanSpec = diff --git a/test/PDBSpec.hs b/test/PDBSpec.hs --- a/test/PDBSpec.hs +++ b/test/PDBSpec.hs @@ -27,33 +27,34 @@ rawPDBToModelConversionSingleChainSpec :: SpecWith () rawPDBToModelConversionSingleChainSpec = describe "Cobot Model from raw single chain PDB" $ do - modelFromPDB <- runIO $ firstPDBModel "test/PDB/1PPE_I.pdb" - modelFromMae <- runIO $ firstMaeModel "test/PDB/1PPE_I.mae" - let (pdbBondCount, _, pdbChainCount, pdbAtomCount) = getStats modelFromPDB - let (maeBondCount, _, _, _) = getStats modelFromMae - - it "Should have correct number of atoms" $ pdbAtomCount `shouldBe` 436 - it "Should have correct number of chains" $ pdbChainCount `shouldBe` 1 - it "Should restore bonds correctly" $ pdbBondCount `shouldBe` maeBondCount + let pdbStats = getStats <$> firstPDBModel "test/PDB/1PPE_I.pdb" + let maeStats = getStats <$> firstMaeModel "test/PDB/1PPE_I.mae" + beforeAll ((,) <$> pdbStats <*> maeStats) $ do + it "Should have correct number of atoms" $ \((_, _, _, pdbAtomCount), _) -> + pdbAtomCount `shouldBe` 436 + it "Should have correct number of chains" $ \((_, _, pdbChainCount, _), _) -> + pdbChainCount `shouldBe` 1 + it "Should restore bonds correctly" $ \((pdbBondCount, _, _, _), (maeBondCount, _, _, _)) -> + pdbBondCount `shouldBe` maeBondCount --- tripeptides are not checked as in checkBiggerMolecule +-- tripeptides are not checked as in checkBiggerMolecule -- because there are inconsistencies in atom numbers between pdb and mae bondsRestoringTripeptideSpec :: SpecWith () bondsRestoringTripeptideSpec = describe "Bonds should be restored correctly in tripeptides" $ sequence_ $ checkTripeptide <$> tripeptides where tripeptides :: [String] - tripeptides = ["ALA_3", "ARG_3", "ASN_3", "ASP_3", "CYS_3", "GLN_3", "GLU_3", "GLY_3", "HID_3", "HIE_3", "HIP_3", + tripeptides = ["ALA_3", "ARG_3", "ASN_3", "ASP_3", "CYS_3", "GLN_3", "GLU_3", "GLY_3", "HID_3", "HIE_3", "HIP_3", "ILE_3", "LEU_3", "LYS_3", "MET_3", "PHE_3", "PRO_3", "SER_3", "THR_3", "TRP_3", "TYR_3", "VAL_3"] checkTripeptide :: String -> SpecWith (Arg Expectation) - checkTripeptide tripeptideName = do - modelFromMae <- runIO . firstMaeModel $ "test/PDB/BondsRestoring/" ++ tripeptideName ++ ".mae" - modelFromPDB <- runIO . firstPDBModel $ "test/PDB/BondsRestoring/" ++ tripeptideName ++ ".pdb" + checkTripeptide tripeptideName = it (tripeptideName ++ " equal bond count in Mae and PDB") $ do + modelFromMae <- firstMaeModel $ "test/PDB/BondsRestoring/" ++ tripeptideName ++ ".mae" + modelFromPDB <- firstPDBModel $ "test/PDB/BondsRestoring/" ++ tripeptideName ++ ".pdb" - let pdbBondCount = V.length $ modelBonds modelFromPDB - let maeBondCount = V.length $ modelBonds modelFromMae + let pdbBondCount = V.length $ modelBonds modelFromPDB + let maeBondCount = V.length $ modelBonds modelFromMae - it (tripeptideName ++ " equal bond count in Mae and PDB") $ pdbBondCount `shouldBe` maeBondCount + pdbBondCount `shouldBe` maeBondCount bondsRestoringBiggerMoleculesSpec :: SpecWith () bondsRestoringBiggerMoleculesSpec = describe "Bonds should be restored correctly in bigger molecules" $ do @@ -62,31 +63,36 @@ checkBiggerMolecule "4dn4_ag_b" where checkBiggerMolecule moleculeName = do - modelFromPDB <- runIO . firstPDBModel $ "test/PDB/BondsRestoring/" ++ moleculeName ++ ".pdb" - modelFromMae <- runIO . firstMaeModel $ "test/PDB/BondsRestoring/" ++ moleculeName ++ ".mae" - let (pdbGlobalBondCount, pdbLocalBondCount, _, _) = getStats modelFromPDB - let (maeGlobalBondCount, maeLocalBondCount, _, _) = getStats modelFromMae - - it (moleculeName ++ " equal global bond count in Mae and PDB") $ pdbGlobalBondCount `shouldBe` maeGlobalBondCount - it (moleculeName ++ " equal local bond count in Mae and PDB") $ pdbLocalBondCount `shouldBe` maeLocalBondCount - - it (moleculeName ++ " no dublicate bonds") $ length (doubleBonds modelFromPDB) `shouldBe` 0 - - let _globalBondSetPDB = globalBondSet modelFromPDB - let _globalBondSetMae = globalBondSet modelFromMae - let diffMaePDBGlobal = S.difference _globalBondSetMae _globalBondSetPDB - let diffPDBMaeGlobal = S.difference _globalBondSetPDB _globalBondSetMae + let pdbIO = firstPDBModel $ "test/PDB/BondsRestoring/" ++ moleculeName ++ ".pdb" + let maeIO = firstMaeModel $ "test/PDB/BondsRestoring/" ++ moleculeName ++ ".mae" + beforeAll ((,) <$> pdbIO <*> maeIO) $ do + it (moleculeName ++ " equal global bond count in Mae and PDB") $ \(modelFromPDB, modelFromMae) -> do + let (pdbGlobalBondCount, _, _, _) = getStats modelFromPDB + let (maeGlobalBondCount, _, _, _) = getStats modelFromMae + pdbGlobalBondCount `shouldBe` maeGlobalBondCount + it (moleculeName ++ " equal local bond count in Mae and PDB") $ \(modelFromPDB, modelFromMae) -> do + let (_, pdbLocalBondCount, _, _) = getStats modelFromPDB + let (_, maeLocalBondCount, _, _) = getStats modelFromMae + pdbLocalBondCount `shouldBe` maeLocalBondCount - it (moleculeName ++ " difference in Mae and PDB global bond sets") $ S.size diffMaePDBGlobal `shouldBe` 0 - it (moleculeName ++ " difference in PDB and Mae global bond sets") $ S.size diffPDBMaeGlobal `shouldBe` 0 + it (moleculeName ++ " no dublicate bonds") $ \(modelFromPDB, _) -> + length (doubleBonds modelFromPDB) `shouldBe` 0 - let _localBondSetPDB = localBondSet modelFromPDB - let _localBondSetMae = localBondSet modelFromMae - let diffMaePDBLocal = S.difference _localBondSetMae _localBondSetPDB - let diffPDBMaeLocal = S.difference _localBondSetPDB _localBondSetMae + it (moleculeName ++ " difference in Mae and PDB global bond sets") $ \(modelFromPDB, modelFromMae) -> do + let _globalBondSetPDB = globalBondSet modelFromPDB + let _globalBondSetMae = globalBondSet modelFromMae + let diffPDBMaeGlobal = S.difference _globalBondSetPDB _globalBondSetMae + let diffMaePDBGlobal = S.difference _globalBondSetMae _globalBondSetPDB + S.size diffPDBMaeGlobal `shouldBe` 0 + S.size diffMaePDBGlobal `shouldBe` 0 - it (moleculeName ++ " difference in Mae and PDB local bond sets") $ S.size diffMaePDBLocal `shouldBe` 0 - it (moleculeName ++ " difference in PDB and Mae local bond sets") $ S.size diffPDBMaeLocal `shouldBe` 0 + it (moleculeName ++ " difference in Mae and PDB local bond sets") $ \(modelFromPDB, modelFromMae) -> do + let _localBondSetPDB = localBondSet modelFromPDB + let _localBondSetMae = localBondSet modelFromMae + let diffMaePDBLocal = S.difference _localBondSetMae _localBondSetPDB + let diffPDBMaeLocal = S.difference _localBondSetPDB _localBondSetMae + S.size diffMaePDBLocal `shouldBe` 0 + S.size diffPDBMaeLocal `shouldBe` 0 localBondSet :: Model -> Set (Text, Int, Int, Int) -- (ChainID, ResidueNumber, LocalFrom, LocalTo) localBondSet Model{..} = S.fromList $ do @@ -127,7 +133,7 @@ getStats :: Model -> (Int, Int, Int, Int) getStats model = (globalBondCount, localBondCount, chainCount, atomCount) - where + where globalBondCount = V.length $ modelBonds model localBondCount = V.length . V.concatMap (V.concatMap resBonds . chainResidues) $ modelChains model chainCount = V.length $ modelChains model diff --git a/test/StructureSpec.hs b/test/StructureSpec.hs --- a/test/StructureSpec.hs +++ b/test/StructureSpec.hs @@ -23,12 +23,12 @@ structureSpec :: Spec structureSpec = describe "Structure spec." $ do - m <- runIO $ V.toList . modelsOf <$> fromFile "test/MAE/Capri.mae" >>= \[x] -> pure x - - it "atoms filtering works correctly. only N, CA, C" $ checkFiltering m $ (`elem` ["N", "CA", "C"]) . atomName - it "atoms filtering works correctly. only CA" $ checkFiltering m $ (== "CA") . atomName - it "atoms filtering works correctly. no atoms" $ checkFiltering m $ const False - it "atoms filtering works correctly. all atoms" $ checkFiltering m $ const True + let mIO = V.toList . modelsOf <$> fromFile "test/MAE/Capri.mae" >>= \[x] -> pure x + beforeAll mIO $ do + it "atoms filtering works correctly. only N, CA, C" $ \m -> checkFiltering m $ (`elem` ["N", "CA", "C"]) . atomName + it "atoms filtering works correctly. only CA" $ \m -> checkFiltering m $ (== "CA") . atomName + it "atoms filtering works correctly. no atoms" $ \m -> checkFiltering m $ const False + it "atoms filtering works correctly. all atoms" $ \m -> checkFiltering m $ const True where checkFiltering :: Model -> (Atom -> Bool) -> Expectation checkFiltering m p = do