diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -2,6 +2,17 @@ ## [Unreleased] +## [0.1.3.0] - 2020-03-27 +### Added +- Residue index in `Structure`. +- Atom input index in `Structure`. +- Bond restoring for PDB. +- Tests for PDB -> Model conversion. +### Changed +- GlobalID now 0-based in mae, PDB, and MMTF. +### Fixed +- A lot of things. + ## [0.1.2.10] - 2020-03-27 ### Added - Lenses for `Structure`. 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: 9e1717465825e1ea9fa611d37fb048cda60ddcea37b8afd16ba77a858416e927 +-- hash: a5a2fbb335e61a22b55a3d74cd08557730ea19a8b93005c693a4576dad82e498 name: cobot-io -version: 0.1.2.10 +version: 0.1.3.0 synopsis: Biological data file formats and IO description: Please see the README on GitHub at category: Bio @@ -50,6 +50,8 @@ Bio.MMTF.MessagePack Bio.MMTF.Type Bio.PDB + Bio.PDB.BondRestoring + Bio.PDB.Functions Bio.PDB.Parser Bio.PDB.Reader Bio.PDB.Type @@ -103,6 +105,7 @@ MAEParserSpec MAESpec MMTFSpec + PDBParserSpec PDBSpec SequenceSpec StructureSpec diff --git a/src/Bio/MAE.hs b/src/Bio/MAE.hs --- a/src/Bio/MAE.hs +++ b/src/Bio/MAE.hs @@ -1,5 +1,5 @@ -{-# LANGUAGE AllowAmbiguousTypes #-} -{-# LANGUAGE CPP #-} +{-# LANGUAGE AllowAmbiguousTypes #-} +{-# LANGUAGE CPP #-} {-# OPTIONS_GHC -fno-warn-orphans #-} module Bio.MAE @@ -8,13 +8,15 @@ , Table (..) , fromFile , fromText + , modelsFromMaeText + , modelsFromMaeFile , maeP ) where import Bio.MAE.Parser import Bio.MAE.Type (Block (..), FromMaeValue (..), Mae (..), MaeValue (..), Table (..)) -import Bio.Structure (Atom (..), Bond (..), Chain (..), +import Bio.Structure (Atom (..), Bond (..), Chain (..), Model (..), GlobalID (..), LocalID (..), Model (..), Residue (..), SecondaryStructure (..), @@ -50,6 +52,12 @@ fromText :: Text -> Either Text Mae fromText = first T.pack . parseOnly maeP +modelsFromMaeFile :: (MonadIO m) => FilePath -> m (Either Text (Vector Model)) +modelsFromMaeFile = liftIO . fmap modelsFromMaeText . TIO.readFile + +modelsFromMaeText :: Text -> Either Text (Vector Model) +modelsFromMaeText maeText = modelsOf <$> fromText maeText + instance StructureModels Mae where modelsOf Mae{..} = V.fromList $ fmap blockToModel blocks where @@ -117,18 +125,23 @@ groupedByResidues = toGroupsOn by group residues = V.fromList $ fmap groupToResidue groupedByResidues - by :: Int -> (Int, Text) - by i = (unsafeGetFromContents "i_m_residue_number" i, getFromContents defaultChainName "s_m_insertion_code" i) + by :: Int -> (Int, Char) + by i = (unsafeGetFromContents "i_m_residue_number" i, getFromContents defaultInsertionCode "s_m_insertion_code" i) defaultChainName :: Text defaultChainName = "A" + defaultInsertionCode :: Char + defaultInsertionCode = ' ' + 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 + groupToResidue group@(h : _) = Residue name residueNumber insertionCode atoms (V.fromList localBonds) secondary chemCompType where - name = stripQuotes $ unsafeGetFromContents "s_m_pdb_residue_name" h - atoms = V.fromList $ fmap indexToAtom group + name = stripQuotes $ unsafeGetFromContents "s_m_pdb_residue_name" h + residueNumber = unsafeGetFromContents "i_m_residue_number" h + insertionCode = unsafeGetFromContents "s_m_insertion_code" h + atoms = V.fromList $ fmap indexToAtom group localInds = [0 .. length group - 1] globalToLocal = M.fromList $ zip group localInds @@ -146,6 +159,7 @@ indexToAtom :: Int -> Atom indexToAtom i = Atom (GlobalID i) + (i + 1) (stripQuotes $ getFromContentsI "s_m_pdb_atom_name") (elIndToElement M.! getFromContentsI "i_m_atomic_number") coords diff --git a/src/Bio/MAE/Type.hs b/src/Bio/MAE/Type.hs --- a/src/Bio/MAE/Type.hs +++ b/src/Bio/MAE/Type.hs @@ -9,6 +9,7 @@ import Data.Map.Strict (Map) import Data.Maybe (fromJust) import Data.Text (Text) +import qualified Data.Text as T (head, null, dropAround) data Mae = Mae { version :: Text , blocks :: [Block] @@ -58,3 +59,11 @@ fromMaeValue :: MaeValue -> Maybe Text fromMaeValue (StringMaeValue t) = Just t fromMaeValue _ = Nothing + +instance FromMaeValue Char where + fromMaeValue :: MaeValue -> Maybe Char + fromMaeValue (StringMaeValue t) = Just $ if T.null t then ' ' else T.head $ stripQuotes t + fromMaeValue _ = Nothing + +stripQuotes :: Text -> Text +stripQuotes = T.dropAround (== '"') diff --git a/src/Bio/MMTF.hs b/src/Bio/MMTF.hs --- a/src/Bio/MMTF.hs +++ b/src/Bio/MMTF.hs @@ -41,6 +41,7 @@ decode (getResponseBody resp) instance StructureModels MMTF where + -- TODO: add global bonds modelsOf m = l2v (flip Model empty . l2v <$> zipWith (zipWith Chain) chainNames chainResis) where chainsCnts = fromIntegral <$> toList (chainsPerModel (model m)) @@ -70,7 +71,8 @@ in (end, mkAtom <$> zip4 cl nl el ics) mkResidue :: (GroupType, SecondaryStructure, [Atom]) -> Residue - mkResidue (gt, ss, atoms') = Residue (gtGroupName gt) (l2v atoms') + -- TODO: support residue number here + mkResidue (gt, ss, atoms') = Residue (gtGroupName gt) (-1) ' ' (l2v atoms') (mkBonds (gtBondAtomList gt) (gtBondOrderList gt)) ss (gtChemCompType gt) @@ -87,13 +89,14 @@ z = zCoordList (atom m) o = occupancyList (atom m) b = bFactorList (atom m) - in Atom (GlobalID $ fromIntegral (i ! idx)) - n - e - (V3 (x ! idx) (y ! idx) (z ! idx)) - fc - (b ! idx) - (o ! idx) + in Atom (GlobalID idx) + (fromIntegral $ i ! idx) + n + e + (V3 (x ! idx) (y ! idx) (z ! idx)) + fc + (b ! idx) + (o ! idx) cutter :: [Int] -> [a] -> [[a]] cutter [] [] = [] diff --git a/src/Bio/PDB.hs b/src/Bio/PDB.hs --- a/src/Bio/PDB.hs +++ b/src/Bio/PDB.hs @@ -1,61 +1,82 @@ {-# OPTIONS_GHC -fno-warn-orphans #-} module Bio.PDB - ( + ( modelsFromPDBText + , modelsFromPDBFile ) where -import qualified Bio.PDB.Type as PDB +import qualified Bio.PDB.Type as PDB +import Bio.PDB.Reader (fromTextPDB, PDBWarnings) +import Bio.PDB.BondRestoring (restoreModelGlobalBonds, restoreChainLocalBonds, residueID) +import Bio.PDB.Functions (groupChainByResidue) import Bio.Structure -import Control.Arrow ((&&&)) -import Data.Coerce (coerce) -import Data.Foldable (Foldable (..)) -import Data.Text as T (Text, singleton, unpack) -import qualified Data.Vector as V -import Linear.V3 (V3 (..)) +import Control.Arrow ((&&&)) +import Control.Monad.IO.Class (MonadIO, liftIO) +import Data.Text as T (Text, singleton, unpack, strip) +import Data.Text.IO as TIO (readFile) +import Data.Map (Map) +import qualified Data.Map as M ((!), fromList) +import qualified Data.Vector as V +import Data.List (sort) +import Data.Maybe (fromMaybe) + +import Text.Read (readMaybe) + +import Linear.V3 (V3 (..)) + instance StructureModels PDB.PDB where modelsOf PDB.PDB {..} = fmap mkModel models where mkModel :: PDB.Model -> Model - mkModel = flip Model V.empty . fmap mkChain - - mkChain :: PDB.Chain -> Chain - mkChain = uncurry Chain . (mkChainName &&& mkChainResidues) + mkModel model = Model (fmap mkChain model) (restoreModelGlobalBonds atomSerialToNilBasedIndex model) + where + atomSerialToNilBasedIndex :: Map Int Int + atomSerialToNilBasedIndex = M.fromList $ allModelAtomSerials `zip` [0..] - mkChainName :: PDB.Chain -> Text - mkChainName = T.singleton . PDB.atomChainID . safeFirstAtom + allModelAtomSerials :: [Int] + allModelAtomSerials = sort . V.toList . fmap PDB.atomSerial . V.concat $ V.toList model - mkChainResidues :: PDB.Chain -> V.Vector Residue - mkChainResidues = V.fromList . fmap mkResidue . flip groupByResidue [] . pure . toList + mkChain :: PDB.Chain -> Chain + mkChain = uncurry Chain . (mkChainName &&& mkChainResidues) - -- can be rewritten with sortOn and groupBy - groupByResidue :: [[PDB.Atom]] -> [PDB.Atom] -> [[PDB.Atom]] - groupByResidue res [] = res - groupByResidue [] (x : xs) = groupByResidue [[x]] xs - groupByResidue res@(lastList : resultTail) (x : xs) - | (PDB.atomResSeq x, PDB.atomICode x) == (PDB.atomResSeq (head lastList), PDB.atomICode (head lastList)) - = groupByResidue ((x : lastList) : resultTail) xs - | otherwise = groupByResidue ([x] : res) xs + mkChainName :: PDB.Chain -> Text + mkChainName = T.singleton . PDB.atomChainID . safeFirstAtom - safeFirstAtom :: V.Vector PDB.Atom -> PDB.Atom - safeFirstAtom arr | V.length arr > 0 = arr V.! 0 - | otherwise = error "Could not pick first atom" + mkChainResidues :: PDB.Chain -> V.Vector Residue + mkChainResidues chain = V.fromList . fmap (mkResidue (restoreChainLocalBonds chain)) $ groupChainByResidue chain + safeFirstAtom :: V.Vector PDB.Atom -> PDB.Atom + safeFirstAtom arr | V.length arr > 0 = arr V.! 0 + | otherwise = error "Could not pick first atom" + + mkResidue :: Map Text (V.Vector (Bond LocalID)) -> [PDB.Atom] -> Residue + mkResidue _ [] = error "Cound not make residue from empty list" + mkResidue localBondsMap atoms' = Residue (T.strip $ PDB.atomResName firstResidueAtom) + (PDB.atomResSeq firstResidueAtom) + (PDB.atomICode firstResidueAtom) + (V.fromList $ mkAtom <$> atoms') + (localBondsMap M.! residueID firstResidueAtom) + Undefined -- now we do not read secondary structure + "" -- chemical component type?! + where + firstResidueAtom = head atoms' - mkResidue :: [PDB.Atom] -> Residue - mkResidue [] = error "Cound not make residue from empty list" - mkResidue atoms' = Residue (PDB.atomResName . head $ atoms') - (V.fromList $ mkAtom <$> atoms') - V.empty -- now we do not read bonds - Undefined -- now we do not read secondary structure - "" -- chemical component type?! + mkAtom :: PDB.Atom -> Atom + mkAtom PDB.Atom{..} = Atom (GlobalID $ atomSerialToNilBasedIndex M.! atomSerial) + atomSerial + (T.strip atomName) + atomElement + (V3 atomX atomY atomZ) + (fromMaybe 0 . readMaybe $ T.unpack atomCharge) + atomTempFactor + atomOccupancy +modelsFromPDBFile :: (MonadIO m) => FilePath -> m (Either Text ([PDBWarnings], V.Vector Model)) +modelsFromPDBFile = liftIO . fmap modelsFromPDBText . TIO.readFile - mkAtom :: PDB.Atom -> Atom - mkAtom PDB.Atom{..} = Atom (coerce atomSerial) - atomName - atomElement - (V3 atomX atomY atomZ) - (read $ T.unpack atomCharge) - atomTempFactor - atomOccupancy +modelsFromPDBText :: Text -> Either Text ([PDBWarnings], V.Vector Model) +modelsFromPDBText pdbText = do + (warnings, parsedPDB) <- fromTextPDB pdbText + let models = modelsOf parsedPDB + pure (warnings, models) diff --git a/src/Bio/PDB/BondRestoring.hs b/src/Bio/PDB/BondRestoring.hs new file mode 100644 --- /dev/null +++ b/src/Bio/PDB/BondRestoring.hs @@ -0,0 +1,196 @@ +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE TupleSections #-} + +module Bio.PDB.BondRestoring + ( restoreModelGlobalBonds + , restoreModelLocalBonds + , restoreChainLocalBonds + , residueID + ) where + +import qualified Bio.PDB.Type as PDB (Atom(..)) +import Bio.PDB.Functions (groupChainByResidue) +import Bio.Structure (Bond (..), GlobalID (..), LocalID (..)) + +import Data.Vector (Vector) +import qualified Data.Vector as V (fromList, toList) +import Data.List (find, sort) +import Data.Text (Text) +import qualified Data.Text as T (strip, pack, unpack) +import Data.Map.Strict (Map, (!?), (!)) +import qualified Data.Map.Strict as M (fromList) +import Data.Maybe (catMaybes) + +import Linear.Metric (distance) +import Linear.V3 (V3(..)) + +import Control.Monad (guard) + +residueID :: PDB.Atom -> Text +residueID PDB.Atom{..} = T.pack (show atomChainID) <> T.pack (show atomResSeq) <> T.pack (show atomICode) + +restoreModelLocalBonds :: Vector (Vector PDB.Atom) -> Map Text (Vector (Bond LocalID)) +restoreModelLocalBonds = M.fromList . concatMap restoreChainLocalBonds' + +restoreChainLocalBonds :: Vector PDB.Atom -> Map Text (Vector (Bond LocalID)) +restoreChainLocalBonds = M.fromList . restoreChainLocalBonds' + +restoreChainLocalBonds' :: Vector PDB.Atom -> [(Text, Vector (Bond LocalID))] +restoreChainLocalBonds' chainAtoms = residueIDToLocalBonds + where + residueIDToLocalBonds :: [(Text, Vector (Bond LocalID))] + residueIDToLocalBonds = do + (residueAtoms, residueBonds) <- zip chainAtomsGroupedByResidue intraResidueGlobalBonds + let localBonds = V.fromList $ convertGlobalsToLocals residueAtoms residueBonds + let _residueID = residueID $ head residueAtoms + pure (_residueID, localBonds) + + intraResidueGlobalBonds :: [[Bond GlobalID]] + intraResidueGlobalBonds = fmap restoreIntraResidueBonds chainAtomsGroupedByResidue + + chainAtomsGroupedByResidue :: [[PDB.Atom]] + chainAtomsGroupedByResidue = groupChainByResidue chainAtoms + + convertGlobalsToLocals :: [PDB.Atom] -> [Bond GlobalID] -> [Bond LocalID] + convertGlobalsToLocals residueAtoms = map convertGlobalToLocal + where + convertGlobalToLocal :: Bond GlobalID -> Bond LocalID + convertGlobalToLocal (Bond (GlobalID from) (GlobalID to) order) = + Bond (LocalID $ globalToLocalIdxMap ! from) (LocalID $ globalToLocalIdxMap ! to) order + + globalToLocalIdxMap :: Map Int Int + globalToLocalIdxMap = M.fromList $ zip sortedGlobalIndices [0..] + + sortedGlobalIndices :: [Int] + sortedGlobalIndices = map PDB.atomSerial $ sort residueAtoms + + +restoreModelGlobalBonds :: Map Int Int -> Vector (Vector PDB.Atom) -> Vector (Bond GlobalID) +restoreModelGlobalBonds atomSerialToNilBasedIndex chains = convertGlobalIDs atomSerialToNilBasedIndex . V.fromList $ _intraResidueBonds ++ peptideBonds ++ disulfideBonds + where + convertGlobalIDs :: Map Int Int -> Vector (Bond GlobalID) -> Vector (Bond GlobalID) + convertGlobalIDs mapping = reindexBonds (\(GlobalID v) -> GlobalID $ mapping ! v) + + reindexBonds :: (a -> a) -> Vector (Bond a) -> Vector (Bond a) + reindexBonds convertID = fmap (\(Bond from to order) -> Bond (convertID from) (convertID to) order) + + chainAtomsGroupedByResidue :: Vector [[PDB.Atom]] + chainAtomsGroupedByResidue = fmap groupChainByResidue chains + + _intraResidueBonds :: [Bond GlobalID] + _intraResidueBonds = concatMap restoreChainIntraResidueBonds chainAtomsGroupedByResidue + + peptideBonds :: [Bond GlobalID] + peptideBonds = concatMap restoreChainPeptideBonds chainAtomsGroupedByResidue + + disulfideBonds :: [Bond GlobalID] + disulfideBonds = restoreDisulfideBonds . concat $ V.toList chainAtomsGroupedByResidue + +restoreDisulfideBonds :: [[PDB.Atom]] -> [Bond GlobalID] +restoreDisulfideBonds atomsGroupedByResidue = do + atom1 <- cystineSulfur + atom2 <- cystineSulfur + guard (PDB.atomSerial atom1 < PDB.atomSerial atom2) + guard $ distance (coords atom1) (coords atom2) < sulfidicBondMaxLength + pure $ Bond (GlobalID $ PDB.atomSerial atom1) (GlobalID $ PDB.atomSerial atom2) 1 + where + cystineSulfur :: [PDB.Atom] + cystineSulfur = filter (("SG" ==) . T.strip . PDB.atomName) $ concat cystines + cystines :: [[PDB.Atom]] + cystines = filter cystinePredicate atomsGroupedByResidue + cystinePredicate :: [PDB.Atom] -> Bool + cystinePredicate residue = any (("SG" ==) . T.strip . PDB.atomName) residue && all (("HG" /=) . T.strip . PDB.atomName) residue + coords :: PDB.Atom -> V3 Float + coords PDB.Atom{..} = V3 atomX atomY atomZ + +sulfidicBondMaxLength :: Float +sulfidicBondMaxLength = 2.56 + +restoreChainPeptideBonds :: [[PDB.Atom]] -> [Bond GlobalID] +restoreChainPeptideBonds atomsGroupedByResidue = restoreChainPeptideBonds' atomsGroupedByResidue [] + where + restoreChainPeptideBonds' :: [[PDB.Atom]] -> [Bond GlobalID] -> [Bond GlobalID] + restoreChainPeptideBonds' [] acc = acc + restoreChainPeptideBonds' [_] acc = acc + restoreChainPeptideBonds' (residue1:residue2:residues) acc = + restoreChainPeptideBonds' (residue2:residues) (constructBond residue1 residue2 : acc) + + constructBond :: [PDB.Atom] -> [PDB.Atom] -> Bond GlobalID + constructBond residue1 residue2 = Bond (GlobalID $ getAtomIndex residue1 "C") (GlobalID $ getAtomIndex residue2 "N") 1 + + getAtomIndex :: [PDB.Atom] -> Text -> Int + getAtomIndex atoms atomNameToFind = case find ((atomNameToFind ==) . T.strip . PDB.atomName) atoms of + Just PDB.Atom{..} -> atomSerial + Nothing -> error ("Atom with name " ++ T.unpack atomNameToFind ++ " wasn't found in residue " ++ residueId atoms ++ ", chain: " ++ chainId atoms) + residueId :: [PDB.Atom] -> String + residueId [] = error "cobot-io: it's impossible to form a residue ID on a residue with no atoms" + residueId (PDB.Atom{..}:_) = T.unpack atomResName ++ show atomResSeq ++ show atomICode + + chainId :: [PDB.Atom] -> String + chainId [] = error "cobot-io: it's impossible to get a chain ID on a chain with no atoms" + chainId (PDB.Atom{..}:_) = show atomChainID + + +restoreChainIntraResidueBonds :: [[PDB.Atom]] -> [Bond GlobalID] +restoreChainIntraResidueBonds = concatMap restoreIntraResidueBonds + +restoreIntraResidueBonds :: [PDB.Atom] -> [Bond GlobalID] +restoreIntraResidueBonds residueAtoms = catMaybes $ constructBond <$> residueBonds + where + -- TODO: support bond order somehow + constructBond :: (Text, Text) -> Maybe (Bond GlobalID) + constructBond (fromAtomName, toAtomName) = Bond <$> constructGlobalID fromAtomName <*> constructGlobalID toAtomName <*> Just 1 + + constructGlobalID :: Text -> Maybe GlobalID + constructGlobalID atomName = GlobalID <$> atomNameToIndex !? atomName + + atomNameToIndex :: Map Text Int + atomNameToIndex = M.fromList $ (\PDB.Atom{..} -> (T.strip atomName, atomSerial)) <$> residueAtoms + + residueBonds :: [(Text, Text)] + residueBonds = intraResidueBonds . T.strip . PDB.atomResName $ head residueAtoms + +intraResidueBonds :: Text -> [(Text, Text)] +intraResidueBonds "NMA" = [("CA", "N")] +intraResidueBonds "ACE" = [("C", "O"), ("C", "CH3")] +intraResidueBonds residueName = backboneBonds ++ caCbBonds residueName ++ sideChainBonds residueName + +backboneBonds :: [(Text, Text)] +backboneBonds = [("N", "CA"), ("CA", "C"), ("C", "O"), ("N", "H")] ++ [("C","OXT"), ("C","HXT")] ++ bwhMany [("N", ["H1", "H2", "H3"])] + +caCbBonds :: Text -> [(Text, Text)] +caCbBonds aminoacid = case aminoacid of + "GLY" -> bwhMany [("CA", ["HA2", "HA3"])] + _ -> [("CA", "CB"), ("CA", "HA")] + +sideChainBonds :: Text -> [(Text, Text)] +sideChainBonds "ALA" = bwhMany [("CB", ["HB1", "HB2", "HB3"])] +sideChainBonds "ARG" = [("CB", "CG"), ("CG", "CD"), ("CD", "NE"), ("NE", "CZ"), ("CZ", "NH2"), ("CZ", "NH1")] ++ bwhMany[("CB", ["HB3", "HB2"]), ("CG", ["HG3", "HG2"]), ("CD", ["HD3", "HD2"]), ("NE", ["HE"]), ("NH1", ["HH12", "HH11"]), ("NH2", ["HH22", "HH21"])] +sideChainBonds "ASN" = [("CB", "CG"), ("CG", "OD1"), ("CG", "ND2")] ++ bwhMany [("CB", ["HB3", "HB2"]), ("ND2", ["HD22", "HD21"])] +sideChainBonds "ASP" = [("CB", "CG"), ("CG", "OD1"), ("CG", "OD2")] ++ bwhMany [("CB", ["HB3", "HB2"]), ("OD2", ["HD2"])] +sideChainBonds "CYS" = [("CB", "SG")] ++ bwhMany [("CB", ["HB3", "HB2"]), ("SG", ["HG"])] +sideChainBonds "GLN" = [("CB", "CG"), ("CG", "CD"), ("CD", "OE1"), ("CD", "NE2")] ++ bwhMany [("CB", ["HB3", "HB2"]), ("CG", ["HG3", "HG2"]), ("NE2", ["HE22", "HE21"])] +sideChainBonds "GLU" = [("CB", "CG"), ("CG", "CD"), ("CD", "OE1"), ("CD", "OE2")] ++ bwhMany [("CB", ["HB3", "HB2"]), ("CG", ["HG3", "HG2"]), ("OE2", ["HE2"])] +sideChainBonds "GLY" = [] -- nothing +sideChainBonds "HIS" = [("CB", "CG"), ("CG", "ND1"), ("ND1", "CE1"), ("CE1", "NE2"), ("NE2", "CD2"), ("CD2", "CG")] ++ bwhMany [("CB", ["HB3", "HB2"]), ("ND1", ["HD1"]), ("CE1", ["HE1"]), ("NE2", ["HE2"]), ("CD2", ["HD2"])] +sideChainBonds "ILE" = [("CB", "CG1"), ("CB", "CG2"), ("CG1", "CD1")] ++ bwhMany [("CB", ["HB"]), ("CG1", ["HG13", "HG12"]), ("CG2", ["HG21", "HG22", "HG23"]), ("CD1", ["HD11", "HD12", "HD13"])] +sideChainBonds "LEU" = [("CB", "CG"), ("CG", "CD1"), ("CG", "CD2")] ++ bwhMany [("CB", ["HB3", "HB2"]), ("CG", ["HG"]), ("CD1", ["HD11", "HD12", "HD13"]), ("CD2", ["HD21", "HD22", "HD23"])] +sideChainBonds "LYS" = [("CB", "CG"), ("CG", "CD"), ("CD", "CE"), ("CE", "NZ")] ++ bwhMany [("CB", ["HB3", "HB2"]), ("CG", ["HG3", "HG2"]), ("CD", ["HD3", "HD2"]), ("CE", ["HE3", "HE2"]), ("NZ", ["HZ1", "HZ2", "HZ3"])] +sideChainBonds "MET" = [("CB", "CG"), ("CG", "SD"), ("SD", "CE")] ++ bwhMany [("CB", ["HB3", "HB2"]), ("CG", ["HG3", "HG2"]), ("CE", ["HE1", "HE2", "HE3"])] +sideChainBonds "PHE" = [("CB", "CG"), ("CG", "CD1"), ("CD1", "CE1"), ("CE1", "CZ"), ("CZ", "CE2"), ("CE2", "CD2"), ("CD2", "CG")] ++ bwhMany [("CB", ["HB3", "HB2"]), ("CD1", ["HD1"]), ("CE1", ["HE1"]), ("CZ", ["HZ"]), ("CE2", ["HE2"]), ("CD2", ["HD2"])] +sideChainBonds "PRO" = [("CB", "CG"), ("CG", "CD"), ("CD", "N")] ++ bwhMany [("CB", ["HB3", "HB2"]), ("CG", ["HG3", "HG2"]), ("CD", ["HD2", "HD3"])] +sideChainBonds "SER" = [("CB", "OG")] ++ bwhMany [("CB", ["HB3", "HB2"]), ("OG", ["HG"])] +sideChainBonds "THR" = [("CB", "OG1"), ("CB", "CG2")] ++ bwhMany [("CB", ["HB"]), ("OG1", ["HG1"]), ("CG2", ["HG21", "HG22", "HG23"])] +sideChainBonds "TRP" = [("CB", "CG"), ("CG", "CD1"), ("CD1", "NE1"), ("NE1", "CE2"), ("CE2", "CD2"), ("CD2", "CG"), ("CD2", "CE3"), ("CE3", "CZ3"), ("CZ3", "CH2"), ("CH2", "CZ2"), ("CZ2", "CE2")] ++ bwhMany [("CB", ["HB3", "HB2"]), ("CD1", ["HD1"]), ("NE1", ["HE1"]), ("CE3", ["HE3"]), ("CZ3", ["HZ3"]), ("CH2", ["HH2"]), ("CZ2", ["HZ2"])] +sideChainBonds "TYR" = [("CB", "CG"), ("CG", "CD1"), ("CD1", "CE1"), ("CE1", "CZ"), ("CZ", "CE2"), ("CE2", "CD2"), ("CD2", "CG"), ("CZ", "OH")] ++ bwhMany [("CB", ["HB3", "HB2"]), ("CD1", ["HD1"]), ("CE1", ["HE1"]), ("CE2", ["HE2"]), ("CD2", ["HD2"]), ("OH", ["HH"])] +sideChainBonds "VAL" = [("CB", "CG1"), ("CB", "CG2")] ++ bwhMany [("CB", ["HB"]), ("CG1", ["HG11", "HG12", "HG13"]), ("CG2", ["HG21", "HG22", "HG23"])] +sideChainBonds unknownResidue = error . T.unpack $ "cobot-io: we don't know what to do with residue " <> unknownResidue + +bwhMany :: [(Text, [Text])] -> [(Text, Text)] +bwhMany = concatMap bwh + +bwh :: (Text, [Text]) -> [(Text, Text)] +bwh = heavyAtomBondsWithHydrogens + +heavyAtomBondsWithHydrogens :: (Text, [Text]) -> [(Text, Text)] +heavyAtomBondsWithHydrogens (heavyAtomName, hydrogenNames) = (heavyAtomName,) <$> hydrogenNames diff --git a/src/Bio/PDB/Functions.hs b/src/Bio/PDB/Functions.hs new file mode 100644 --- /dev/null +++ b/src/Bio/PDB/Functions.hs @@ -0,0 +1,20 @@ +module Bio.PDB.Functions + ( groupChainByResidue + ) where + +import qualified Bio.PDB.Type as PDB (Atom (..)) +import Data.Map (Map) +import qualified Data.Map as M (fromList, (!)) +import Data.List (groupBy, sortOn) +import Data.Vector (Vector) +import qualified Data.Vector as V (toList) + +groupChainByResidue :: Vector PDB.Atom -> [[PDB.Atom]] +groupChainByResidue = sortOn (sortOnResidue . head) . groupBy atomsFromSameResidue . V.toList + where + atomsFromSameResidue :: PDB.Atom -> PDB.Atom -> Bool + atomsFromSameResidue atom1 atom2 = PDB.atomResSeq atom1 == PDB.atomResSeq atom2 && PDB.atomICode atom1 == PDB.atomICode atom2 + sortOnResidue :: PDB.Atom -> Int + sortOnResidue PDB.Atom{..} = atomSerial * 100 + (insertionCodeSortingCorrections M.! atomICode) + insertionCodeSortingCorrections :: Map Char Int + insertionCodeSortingCorrections = M.fromList $ zip (' ':['A'..'Z']) [0..] diff --git a/src/Bio/PDB/Parser.hs b/src/Bio/PDB/Parser.hs --- a/src/Bio/PDB/Parser.hs +++ b/src/Bio/PDB/Parser.hs @@ -64,8 +64,13 @@ atomP :: Parser CoordLike atomP = let atom = Atom <$> - (string "ATOM " *> -- (1 - 6) # we increased atomSerial field for one symbol - (read <$> count 6 notEndLineChar) <* space) -- (7 - 11) atomSerial + ( + (string "ATOM " *> -- (1 - 5) ATOM -- we extended atomSerial length to the left for one symbol + (read <$> count 6 notEndLineChar) <* space) -- (6 - 11) atomSerial + <|> -- or + (string "HETATM" *> -- (1 - 6) HETATM + (read <$> count 5 notEndLineChar) <* space) -- (7 - 11) atomSerial + ) <*> (T.pack <$> count 4 notEndLineChar) -- (13 - 16) atomName <*> notEndLineChar -- (17) atomAltLoc <*> (T.pack <$> count 3 notEndLineChar) <* space -- (18 - 20) atomResName @@ -84,7 +89,7 @@ coordNotAtomP :: Parser CoordLike coordNotAtomP = do - _ <- string "HETATM" <|> string "TER " <|> string "ANISOU" <|> string "CONECT" + _ <- string "TER " <|> string "ANISOU" <|> string "CONECT" skipWhile $ not . isEndOfLine endOfLine return CoordNotAtomLine diff --git a/src/Bio/PDB/Reader.hs b/src/Bio/PDB/Reader.hs --- a/src/Bio/PDB/Reader.hs +++ b/src/Bio/PDB/Reader.hs @@ -1,6 +1,7 @@ module Bio.PDB.Reader ( fromTextPDB , fromFilePDB + , PDBWarnings(..) ) where import Bio.PDB.Parser (pdbP) @@ -61,13 +62,13 @@ then Right standardizedText else Left "There are trash strings between model strings" -fromFilePDB :: MonadIO m => FilePath -> m (Either Text ([PDBWarnings], Either Text PDB)) -fromFilePDB f = do - content <- liftIO (TIO.readFile f) - let preprocessed = preprocess content - pure $ fmap (first pack . parseOnly pdbP) <$> preprocessed -fromTextPDB :: Text -> Either Text ([PDBWarnings], Either Text PDB) -fromTextPDB text = fmap (first pack . parseOnly pdbP) <$> preprocessed - where - preprocessed = preprocess text +fromFilePDB :: MonadIO m => FilePath -> m (Either Text ([PDBWarnings], PDB)) +fromFilePDB = liftIO . fmap fromTextPDB . TIO.readFile + +fromTextPDB :: Text -> Either Text ([PDBWarnings], PDB) +fromTextPDB text = do + (warnings, preprocessedText) <- preprocess text + pdb <- first pack $ parseOnly pdbP preprocessedText + + pure (warnings, pdb) diff --git a/src/Bio/PDB/Type.hs b/src/Bio/PDB/Type.hs --- a/src/Bio/PDB/Type.hs +++ b/src/Bio/PDB/Type.hs @@ -92,3 +92,6 @@ , atomCharge :: Text -- ^ Charge on the atom. } deriving (Show, Eq, Generic, NFData) + +instance Ord Atom where + a1 <= a2 = atomSerial a1 <= atomSerial a2 diff --git a/src/Bio/Structure.hs b/src/Bio/Structure.hs --- a/src/Bio/Structure.hs +++ b/src/Bio/Structure.hs @@ -44,16 +44,20 @@ -- | Generic atom representation -- -data Atom = Atom { atomId :: GlobalID -- ^ global identifier - , atomName :: Text -- ^ IUPAC atom name - , atomElement :: Text -- ^ atom chemical element - , atomCoords :: V3 Float -- ^ 3D coordinates of atom - , formalCharge :: Int -- ^ Formal charge of atom - , bFactor :: Float -- ^ B-factor of atom - , occupancy :: Float -- ^ the amount of each conformation that is observed in the crystal +data Atom = Atom { atomId :: GlobalID -- ^ global identifier, 0-based + , atomInputIndex :: Int -- ^ atom index from input file + , atomName :: Text -- ^ IUPAC atom name + , atomElement :: Text -- ^ atom chemical element + , atomCoords :: V3 Float -- ^ 3D coordinates of atom + , formalCharge :: Int -- ^ Formal charge of atom + , bFactor :: Float -- ^ B-factor of atom + , occupancy :: Float -- ^ the amount of each conformation that is observed in the crystal } deriving (Show, Eq, Generic) +instance Ord Atom where + a1 <= a2 = atomId a1 <= atomId a2 + instance NFData Atom -- | Generic chemical bond @@ -76,11 +80,13 @@ -- | A set of atoms, organized to a residues -- -data Residue = Residue { resName :: Text -- ^ residue name - , resAtoms :: Vector Atom -- ^ a set of residue atoms - , resBonds :: Vector (Bond LocalID) -- ^ a set of residue bonds with local identifiers (position in 'resAtoms') - , resSecondary :: SecondaryStructure -- ^ residue secondary structure - , resChemCompType :: Text -- ^ chemical component type +data Residue = Residue { resName :: Text -- ^ residue name + , resNumber :: Int -- ^ residue number + , resInsertionCode :: Char -- ^ residue insertion code + , resAtoms :: Vector Atom -- ^ a set of residue atoms + , resBonds :: Vector (Bond LocalID) -- ^ a set of residue bonds with local identifiers (position in 'resAtoms') + , resSecondary :: SecondaryStructure -- ^ residue secondary structure + , resChemCompType :: Text -- ^ chemical component type } deriving (Show, Eq, Generic, NFData) diff --git a/test/MAESpec.hs b/test/MAESpec.hs --- a/test/MAESpec.hs +++ b/test/MAESpec.hs @@ -5,86 +5,107 @@ module MAESpec where -import Bio.MAE (fromFile) +import Bio.MAE (modelsFromMaeFile) import Bio.Structure (Atom (..), Bond (..), Chain (..), GlobalID (..), - LocalID (..), Model (..), Residue (..), - StructureModels (..)) + LocalID (..), Model (..), Residue (..)) + +import Control.Monad.IO.Class (MonadIO) + import Data.Set (Set) import qualified Data.Set as S (fromList) import Data.Vector (Vector) -import qualified Data.Vector as V (fromList, toList, (!)) +import qualified Data.Vector as V (fromList, toList, (!), head) +import Data.Either (fromRight) 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 +maeSpec = describe "Mae spec." $ + beforeAll (firstMaeModel "test/MAE/small.mae") $ do + let firstChainResidues Model{..} = (chainResidues $ modelChains V.! 0) + let secondChainResidues Model{..} = (chainResidues $ modelChains V.! 1) - 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 "two chains" $ \Model{..} -> length modelChains `shouldBe` 2 + + it "residue numbers" $ \m@Model{..} -> do + fmap resName (firstChainResidues m) `shouldBe` V.fromList ["ACE", "ASP", "ILE", "LYS"] + fmap resName (secondChainResidues m) `shouldBe` V.fromList ["GLU", "LEU", "VAL", "ARG", "PRO", "GLY", "ALA", "LEU", "VAL"] + + it "residue names" $ \m@Model{..} -> do + fmap resNumber (firstChainResidues m) `shouldBe` V.fromList [0, 1, 2, 3] + fmap resNumber (secondChainResidues m) `shouldBe` V.fromList [10, 11, 12, 13, 13, 15, 16, 17, 18] -- 13 is doubled because the second 13 has 'A' insertion code - it "atoms count" $ do - sum (fmap (length . resAtoms) $ chainResidues $ modelChains V.! 0) `shouldBe` 62 - sum (fmap (length . resAtoms) $ chainResidues $ modelChains V.! 1) `shouldBe` 140 + it "residue insertion codes" $ \m@Model{..} -> do + fmap resInsertionCode (firstChainResidues m) `shouldBe` V.fromList [' ', ' ', ' ', ' '] + fmap resInsertionCode (secondChainResidues m) `shouldBe` V.fromList [' ', ' ', ' ', ' ', 'A', ' ', ' ', ' ', ' '] - 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 "atoms count" $ \m@Model{..} -> do + sum (fmap (length . resAtoms) (firstChainResidues m)) `shouldBe` 62 + sum (fmap (length . resAtoms) (secondChainResidues m)) `shouldBe` 140 - it "global bonds" $ do - modelBonds `shouldBe` toBond GlobalID <$> V.fromList allBonds + 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) + ] - let residue1 = (chainResidues $ modelChains V.! 0) V.! 3 - let residue2 = (chainResidues $ modelChains V.! 1) V.! 0 + it "global bonds" $ \Model{..} -> + modelBonds `shouldBe` toBond GlobalID <$> V.fromList allBonds - it "atoms in residues" $ do - let atoms1 = resAtoms residue1 - let atoms2 = resAtoms residue2 + let residue1 m = firstChainResidues m V.! 3 + let residue2 m = secondChainResidues m V.! 0 - fmap atomId atoms1 `shouldBe` GlobalID <$> V.fromList ([37 .. 58] <> [199 .. 201]) - fmap atomId atoms2 `shouldBe` GlobalID <$> V.fromList [59 .. 73] + it "atoms in residues" $ \m@Model{..} -> do + let atoms1 = resAtoms $ residue1 m + let atoms2 = resAtoms $ residue2 m - 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 atomId atoms1 `shouldBe` GlobalID <$> V.fromList ([37 .. 58] <> [199 .. 201]) + fmap atomId atoms2 `shouldBe` GlobalID <$> V.fromList [59 .. 73] - 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 - ] + 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" + ] - 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 + 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 + ] - toSet :: Ord a => Vector a -> Set a - toSet = S.fromList . V.toList + it "bonds in residues" $ \m@Model{..} -> do + toSet (resBonds $ residue1 m) `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 m) `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 + +firstMaeModel :: (MonadIO m) => FilePath -> m Model +firstMaeModel filepath = do + eitherMae <- modelsFromMaeFile filepath + -- `evaluate . force` fails for some reason + pure . V.head $ fromRight undefined eitherMae diff --git a/test/MMTFSpec.hs b/test/MMTFSpec.hs --- a/test/MMTFSpec.hs +++ b/test/MMTFSpec.hs @@ -2,8 +2,9 @@ import Bio.MMTF import Bio.MMTF.Decode.Codec -import Data.Int (Int8) -import Data.Vector ((!)) +import Data.Int (Int8) +import Data.Vector ((!)) +import Data.ByteString.Lazy as BS (readFile) import Test.Hspec mmtfCodecSpec :: Spec @@ -26,11 +27,10 @@ mmtfParserSpec = describe "MMTF parser" $ it "should parse 1FSD" $ do - m <- fetch "1FSD" + m <- BS.readFile "test/MMTF/1FSD.dms" >>= decode (structureId . structure) m `shouldBe` "1FSD" (numModels . structure) m `shouldBe` 41 (length . bFactorList . atom) m `shouldBe` 20664 ((! 0) . experimentalMethods . structure) m `shouldBe` "SOLUTION NMR" ((! 0) . xCoordList . atom) m `shouldBe` (-12.847) ((! 20663) . xCoordList . atom) m `shouldBe` 5.672 - diff --git a/test/PDBParserSpec.hs b/test/PDBParserSpec.hs new file mode 100644 --- /dev/null +++ b/test/PDBParserSpec.hs @@ -0,0 +1,446 @@ +{-# LANGUAGE OverloadedStrings #-} + +module PDBParserSpec where + +import Bio.PDB.Reader (fromTextPDB) +import Bio.PDB.Type (Atom (..), FieldType (..), PDB (..)) +import qualified Data.Map.Strict (empty, fromList, singleton) +import Data.Text as T (Text, intercalate, length, lines, pack, + replicate, take) +import qualified Data.Vector as V (empty, fromList, singleton) +import Test.Hspec + +oneModelSpecP :: Spec +oneModelSpecP = describe "One model." $ + it "correctly parses pdb with only one model without strings \"MODEL\" & \"ENDMDL\"" $ do + let mt = fromTextPDB . lenghtenLines $ T.pack ( "HEADER header\n" ++ + "TITLE STRUCTURE OF THE TRANSFORMED MONOCLINIC LYSOZYME\n" ++ + "COMPND compnd\n" ++ + "SOURCE source\n" ++ + "KEYWDS keywds\n" ++ + "AUTHOR Masha\n" ++ + "REVDAT revdat\n" ++ + "REMARK 1 REFERENCE 1\n" ++ + "SEQRES seqres\n" ++ + "CRYST1 cryst1\n" ++ + "ORIGX1 origx1 n=1\n" ++ + "SCALE2 sclaen n=2\n" ++ + "ATOM 2032 OXT CYS A 214 -4.546 -29.673 26.796 1.0 143.51 O \n" ++ + "ATOM 2033 H CYS A 214 -6.124 -27.225 26.558 1.00 15.00 H \n" ++ + "TER 2034 CYS A 214 \n" ++ + "ATOM 2035 N GLU B 1 18.637 -61.583 66.852 1.0 118.48 N \n" ++ + "TER 12534 ARG D 474 \n" ++ + "HETATM12535 C1 NAG B 475 5.791 -20.194 -7.051 1.00 34.66 C \n" ++ + "HETATM12538 C4 NAG B 475 6.943 -19.507 -9.597 1.00 25.87 C \n" ++ + "ATOM 2036 CA GLU B 1 19.722 -62.606 66.868 1.00 19.77 C \n" ++ + "ATOM 2036 CA GLU B 1A 19.722 -62.606 66.868 1.00 19.77 C \n" ++ + "CONECT conect\n" ++ + "CONECT conect conect\n" ++ + "MASTER 1 2 3 4 5 6 7 8\n" ++ + "END" + ) + mt `shouldBe` Right ([], oneModelPDB) + +oneModelPDB :: PDB +oneModelPDB = PDB { title = "STRUCTURE OF THE TRANSFORMED MONOCLINIC LYSOZYME" + , remarks = Data.Map.Strict.singleton (Just 1) (V.singleton "REFERENCE 1") + , models = V.singleton $ V.fromList [V.fromList [ + Atom {atomSerial = 2032, atomName = " OXT", atomAltLoc = ' ', atomResName = "CYS", atomChainID = 'A', atomResSeq = 214, atomICode = ' ', + atomX = -4.546, atomY = -29.673, atomZ = 26.796, atomOccupancy = 1.0, atomTempFactor = 143.51, atomElement = " O", atomCharge = " "}, + Atom {atomSerial = 2033, atomName = " H ", atomAltLoc = ' ', atomResName = "CYS", atomChainID = 'A', atomResSeq = 214, atomICode = ' ', + atomX = -6.124, atomY = -27.225, atomZ = 26.558, atomOccupancy = 1.0, atomTempFactor = 15.0, atomElement = " H", atomCharge = " "}], + V.fromList [ + Atom {atomSerial = 2035, atomName = " N ", atomAltLoc = ' ', atomResName = "GLU", atomChainID = 'B', atomResSeq = 1, atomICode = ' ', + atomX = 18.637, atomY = -61.583, atomZ = 66.852, atomOccupancy = 1.0, atomTempFactor = 118.48, atomElement = " N", atomCharge = " "}, + Atom {atomSerial = 12535, atomName = " C1 ", atomAltLoc = ' ', atomResName = "NAG", atomChainID = 'B', atomResSeq = 475, atomICode = ' ', + atomX = 5.791, atomY = -20.194, atomZ = -7.051, atomOccupancy = 1.0, atomTempFactor = 34.66, atomElement = " C", atomCharge = " "}, + Atom {atomSerial = 12538, atomName = " C4 ", atomAltLoc = ' ', atomResName = "NAG", atomChainID = 'B', atomResSeq = 475, atomICode = ' ', + atomX = 6.943, atomY = -19.507, atomZ = -9.597, atomOccupancy = 1.0, atomTempFactor = 25.87, atomElement = " C", atomCharge = " "}, + Atom {atomSerial = 2036, atomName = " CA ", atomAltLoc = ' ', atomResName = "GLU", atomChainID = 'B', atomResSeq = 1, atomICode = ' ', + atomX = 19.722, atomY = -62.606, atomZ = 66.868, atomOccupancy = 1.0, atomTempFactor = 19.77, atomElement = " C", atomCharge = " "}, + Atom {atomSerial = 2036, atomName = " CA ", atomAltLoc = ' ', atomResName = "GLU", atomChainID = 'B', atomResSeq = 1, atomICode = 'A', + atomX = 19.722, atomY = -62.606, atomZ = 66.868, atomOccupancy = 1.0, atomTempFactor = 19.77, atomElement = " C", atomCharge = " "}] + ] + , otherFields = Data.Map.Strict.fromList [(HEADER, V.fromList [" header"]),(COMPND,V.fromList[" compnd"]),(SOURCE,V.fromList[" source"]),(KEYWDS,V.fromList[" keywds"]),(AUTHOR,V.fromList[" Masha"]),(REVDAT,V.fromList[" revdat"]),(SEQRES,V.fromList [" seqres"]), + (CRYST1,V.fromList [" cryst1"]),(ORIGX1,V.fromList[" origx1 n=1"]),(SCALE2,V.fromList [" sclaen n=2"]),(MASTER,V.fromList[" 1 2 3 4 5 6 7 8"])] + } + +manyModelsSpecP :: Spec +manyModelsSpecP = describe "Some models." $ + it "correctly parses pdb with many models - they have strings \"MODEL\" & \"ENDMDL\" and text has disordered strings" $ do + let mt = fromTextPDB . lenghtenLines $ T.pack ( "HEADER header\n" ++ + "TITLE STRUCTURE OF THE TRANSFORMED MONOCLINIC LYSOZYME\n" ++ + "COMPND compnd\n" ++ + "SOURCE source\n" ++ + "REMARK 1 REFERENCE 1\n" ++ + "REMARK 2 Reference 2_1/1\n" ++ + "CRYST1 cryst1\n" ++ + "ORIGX1 origx1 n=1\n" ++ + "SCALE2 sclaen n=2\n" ++ + "MODEL 4 \n" ++ + "ATOM 2033 H CYS A 214 -6.124 -27.225 26.558 1.00 15.00 H \n" ++ + "TER 2034 CYS A 214 \n" ++ + "ANISOU anisou\n" ++ + "ENDMDL \n" ++ + "MODEL 5 \n" ++ + "ATOM 2035 N GLU B 1 18.637-691.583 66.852 1.0 118.48 N \n" ++ + "ATOM 12531 HH12 ARG D 474 45.558 -39.551 -49.936 1.00 15.00 H \n" ++ + "ENDMDL \n" ++ + "MODEL 6 \n" ++ + "ATOM 12532 HH21 ARG D 474 47.457 -38.007 -47.445 1.00 15.00 H \n" ++ + "TER 12534 ARG D 474 \n" ++ + "ATOM 12533 HH22 ARG D 474 47.405 -39.268 -48.629 1.00 15.00 H \n" ++ + "HETATM12535 C1 NAG B 475 5.791 -20.194 -7.051 1.00 34.66 C \n" ++ + "HETATM12538 C4 NAG B 475 6.943 -19.507 -9.597 1.00 25.87 C \n" ++ + "CONECT conect conect\n" ++ + "CONECT conect\n" ++ + "ENDMDL\n" ++ + "SEQRES seqres\n" ++ + "KEYWDS keywds\n" ++ + "EXPDTA expdta\n" ++ + "AUTHOR Masha\n" ++ + "REVDAT revdat\n" ++ + "MASTER 1 2 3 4 5 6 7 8\n" + ) + mt `shouldBe` Right ([], manyModelsPDB) + +manyModelsPDB :: PDB +manyModelsPDB = PDB { title = "STRUCTURE OF THE TRANSFORMED MONOCLINIC LYSOZYME" + , remarks = Data.Map.Strict.fromList [(Just 1, V.singleton "REFERENCE 1"), (Just 2, V.singleton "Reference 2_1/1")] + , models = V.fromList [V.fromList [V.fromList [ + Atom {atomSerial = 2033, atomName = " H ", atomAltLoc = ' ', atomResName = "CYS", atomChainID = 'A', atomResSeq = 214, atomICode = ' ', + atomX = -6.124, atomY = -27.225, atomZ = 26.558, atomOccupancy = 1.0, atomTempFactor = 15.0, atomElement = " H", atomCharge = " "}] + ], + V.fromList [V.fromList [ + Atom {atomSerial = 2035, atomName = " N ", atomAltLoc = ' ', atomResName = "GLU", atomChainID = 'B', atomResSeq = 1, atomICode = ' ', + atomX = 18.637, atomY = -691.583, atomZ = 66.852, atomOccupancy = 1.0, atomTempFactor = 118.48, atomElement = " N", atomCharge = " "}], + V.fromList [ + Atom {atomSerial = 12531, atomName = "HH12", atomAltLoc = ' ', atomResName = "ARG", atomChainID = 'D', atomResSeq = 474, atomICode = ' ', + atomX = 45.558, atomY = -39.551, atomZ = -49.936, atomOccupancy = 1.0, atomTempFactor = 15.0, atomElement = " H", atomCharge = " "}] + ], + V.fromList [V.fromList [ + Atom {atomSerial = 12532, atomName = "HH21", atomAltLoc = ' ', atomResName = "ARG", atomChainID = 'D', atomResSeq = 474, atomICode = ' ', + atomX = 47.457, atomY = -38.007, atomZ = -47.445, atomOccupancy = 1.0, atomTempFactor = 15.0, atomElement = " H", atomCharge = " "}, + Atom {atomSerial = 12533, atomName = "HH22", atomAltLoc = ' ', atomResName = "ARG", atomChainID = 'D', atomResSeq = 474, atomICode = ' ', + atomX = 47.405, atomY = -39.268, atomZ = -48.629, atomOccupancy = 1.0, atomTempFactor = 15.0, atomElement = " H", atomCharge = " "}], + V.fromList [ + Atom {atomSerial = 12535, atomName = " C1 ", atomAltLoc = ' ', atomResName = "NAG", atomChainID = 'B', atomResSeq = 475, atomICode = ' ', + atomX = 5.791, atomY = -20.194, atomZ = -7.051, atomOccupancy = 1.0, atomTempFactor = 34.66, atomElement = " C", atomCharge = " "}, + Atom {atomSerial = 12538, atomName = " C4 ", atomAltLoc = ' ', atomResName = "NAG", atomChainID = 'B', atomResSeq = 475, atomICode = ' ', + atomX = 6.943, atomY = -19.507, atomZ = -9.597, atomOccupancy = 1.0, atomTempFactor = 25.87, atomElement = " C", atomCharge = " "}] + ] + ] + , otherFields = Data.Map.Strict.fromList [(HEADER, V.fromList [" header"]),(COMPND,V.fromList[" compnd"]),(SOURCE,V.fromList[" source"]),(KEYWDS,V.fromList[" keywds"]),(EXPDTA,V.fromList[" expdta"]),(AUTHOR,V.fromList[" Masha"]),(REVDAT,V.fromList[" revdat"]),(SEQRES,V.fromList [" seqres"]), + (CRYST1,V.fromList [" cryst1"]),(ORIGX1,V.fromList[" origx1 n=1"]),(SCALE2,V.fromList [" sclaen n=2"]),(MASTER,V.fromList[" 1 2 3 4 5 6 7 8"])] + } + +noModelsSpecP :: Spec +noModelsSpecP = describe "No models." $ + it "correctly parses pdb without models (no ATOM, TER, HETATM strings)" $ do + let mt = fromTextPDB . lenghtenLines $ T.pack ( "HEADER header\n" ++ + "TITLE STRUCTURE OF THE TRANSFORMED MONOCLINIC LYSOZYME\n" ++ + "COMPND compnd\n" ++ + "SOURCE source\n" ++ + "KEYWDS keywds\n" ++ + "EXPDTA expdta\n" ++ + "AUTHOR Masha\n" ++ + "REVDAT revdat\n" ++ + "REMARK 1 REFERENCE 1\n" ++ + "SEQRES seqres\n" ++ + "CRYST1 cryst1\n" ++ + "ORIGX1 origx1 n=1\n" ++ + "SCALE2 sclaen n=2\n" ++ + "MASTER 1 2 3 4 5 6 7 8\n" + ) + mt `shouldBe` Right ([], noModelsPDB) + +noModelsPDB :: PDB +noModelsPDB = PDB { title = "STRUCTURE OF THE TRANSFORMED MONOCLINIC LYSOZYME" + , remarks = Data.Map.Strict.singleton (Just 1) (V.singleton "REFERENCE 1") + , models = V.empty + , otherFields = Data.Map.Strict.fromList [(HEADER, V.fromList [" header"]),(COMPND,V.fromList[" compnd"]),(SOURCE,V.fromList[" source"]),(KEYWDS,V.fromList[" keywds"]),(EXPDTA,V.fromList[" expdta"]),(AUTHOR,V.fromList[" Masha"]),(REVDAT,V.fromList[" revdat"]),(SEQRES,V.fromList [" seqres"]), + (CRYST1,V.fromList [" cryst1"]),(ORIGX1,V.fromList[" origx1 n=1"]),(SCALE2,V.fromList [" sclaen n=2"]),(MASTER,V.fromList[" 1 2 3 4 5 6 7 8"])] + } + +allFieldsModelSpecP :: Spec +allFieldsModelSpecP = describe "PDB with all strings." $ + it "correctly parses pdb with all types of string" $ do + let mt = fromTextPDB . lenghtenLines $ T.pack ( "HEADER header\n" ++ + "OBSLTE obslte\n" ++ + "TITLE STRUCTURE OF THE TRANSFORMED MONOCLINIC LYSOZYME\n" ++ + "SPLIT split\n" ++ + "CAVEAT caveat\n" ++ + "COMPND compnd\n" ++ + "SOURCE source\n" ++ + "KEYWDS keywds\n" ++ + "EXPDTA expdta1\n" ++ + "EXPDTA expdta2\n" ++ + "NUMMDL nummdl\n" ++ + "MDLTYP mdltyp mdltyp\n" ++ + "AUTHOR Masha\n" ++ + "REVDAT revdat\n" ++ + "SPRSDE sprsde\n" ++ + "JRNL jrnl\n" ++ + "REMARK 1 REFERENCE 1\n" ++ + "DBREF dbref\n" ++ + "DBREF1 dbref1\n" ++ + "DBREF2 dbref2_1\n" ++ + "DBREF2 dbref2_2\n" ++ + "SEQADV seqadv\n" ++ + "SEQRES seqres\n" ++ + "MODRES modres\n" ++ + "HET het\n" ++ + "HETNAM hetnam\n" ++ + "HETSYN hetsyn\n" ++ + "FORMUL of love\n" ++ + "HELIX helix\n" ++ + "SHEET sheet\n" ++ + "SSBOND ssbond\n" ++ + "LINK link\n" ++ + "CISPEP cispep\n" ++ + "SITE site\n" ++ + "CRYST1 cryst1\n" ++ + "ORIGX1 origx1 n=1\n" ++ + "SCALE2 sclaen n=2\n" ++ + "MTRIX3 matrixn n=3\n" ++ + "MODEL 1\n" ++ + "ATOM 12532 HH21 ARG D 474 47.457 -38.007 -47.445 1.00 15.00 H \n" ++ + "TER 12534 ARG D 474 \n" ++ + "HETATM12535 C1 NAG B 475 5.791 -20.194 -7.051 1.00 34.66 C \n" ++ + "CONECT conect\n" ++ + "CONECT conect conect\n" ++ + "ENDMDL\n" ++ + "MASTER 1 2 3 4 5 6 7 8\n" ++ + "END" + ) + mt `shouldBe` Right ([], pdbWithAllFields) + +pdbWithAllFields :: PDB +pdbWithAllFields = PDB { title = "STRUCTURE OF THE TRANSFORMED MONOCLINIC LYSOZYME" + , models = V.fromList [V.fromList [ + V.singleton Atom {atomSerial = 12532, atomName = "HH21", atomAltLoc = ' ', atomResName = "ARG", atomChainID = 'D', atomResSeq = 474, atomICode = ' ', + atomX = 47.457, atomY = -38.007, atomZ = -47.445, atomOccupancy = 1.0, atomTempFactor = 15.0, atomElement = " H", atomCharge = " "}, + V.singleton Atom {atomSerial = 12535, atomName = " C1 ", atomAltLoc = ' ', atomResName = "NAG", atomChainID = 'B', atomResSeq = 475, atomICode = ' ', + atomX = 5.791, atomY = -20.194, atomZ = -7.051, atomOccupancy = 1.0, atomTempFactor = 34.66, atomElement = " C", atomCharge = " "}] + ] + , remarks = Data.Map.Strict.singleton (Just 1) (V.singleton "REFERENCE 1") + , otherFields = Data.Map.Strict.fromList [(HEADER, V.fromList [" header"]),(OBSLTE, V.fromList [" obslte"]), (SPLIT, V.fromList ["split"]), (CAVEAT, V.fromList [" caveat"]), + (COMPND, V.fromList [" compnd"]),(SOURCE, V.fromList [" source"]),(KEYWDS, V.fromList [" keywds"]), + (EXPDTA, V.fromList [" expdta1"," expdta2"]),(NUMMDL, V.fromList [" nummdl"]),(MDLTYP, V.fromList [" mdltyp mdltyp"]), + (AUTHOR, V.fromList [" Masha"]),(REVDAT, V.fromList [" revdat"]),(SPRSDE, V.fromList [" sprsde"]), + (JRNL, V.fromList [" jrnl"]),(DBREF, V.fromList ["dbref"]), (DBREF1, V.fromList [" dbref1"]), (DBREF2, V.fromList [" dbref2_1"," dbref2_2"]), + (SEQADV, V.fromList [" seqadv"]),(SEQRES, V.fromList [" seqres"]),(MODRES, V.fromList [" modres"]), + (HET, V.fromList [" het"]), (HETNAM, V.fromList [" hetnam"]), (HETSYN, V.fromList [" hetsyn"]),(FORMUL, V.fromList [" of love"]),(HELIX, V.fromList ["helix"]), + (SHEET, V.fromList ["sheet"]),(SSBOND, V.fromList [" ssbond"]),(LINK, V.fromList [" link"]), + (CISPEP, V.fromList [" cispep"]),(SITE,V.fromList [" site"]),(CRYST1, V.fromList [" cryst1"]),(MTRIX3, V.fromList [" matrixn n=3"]), + (ORIGX1, V.fromList [" origx1 n=1"]),(SCALE2, V.fromList [" sclaen n=2"]),(MASTER, V.fromList [" 1 2 3 4 5 6 7 8"])] + } + + +emptySpecP :: Spec +emptySpecP = describe "empty PDB." $ + it "correctly parses empty pdb" $ do + let mt = fromTextPDB "" + mt `shouldBe` Right ([], emptyPdb) + +emptyPdb :: PDB +emptyPdb = PDB { title = "" + , models = V.empty + , remarks = Data.Map.Strict.empty + , otherFields = Data.Map.Strict.empty + } + +trashBetweenModelsSpecP :: Spec +trashBetweenModelsSpecP = describe "PDB has trash." $ + it "correctly parses pdb with trash string between models and other field strings" $ do + let mt = fromTextPDB . lenghtenLines $ T.pack ( "trash strings 1\n" ++ + "HEADER header\n" ++ + "TITLE STRUCTURE OF THE TRANSFORMED MONOCLINIC LYSOZYME\n" ++ + "REMARK 1 REFERENCE 1\n" ++ + "trash strings 2\n" ++ + "REMARK 2 Reference 2_1/1\n" ++ + "MODEL 4 \n" ++ + "ATOM 2033 H CYS A 214 -6.124 -27.225 26.558 1.00 15.00 H \n" ++ + "TER 2034 CYS A 214 \n" ++ + "ANISOU anisou\n" ++ + "ENDMDL \n" ++ + "CRYST1 cryst1\n" ++ + "trash strings 3\n" ++ + "SEQRES seqres\n" ++ + "MODEL 5 \n" ++ + "ATOM 2035 N GLU B 1 18.637-691.583 66.852 1.0 118.48 N \n" ++ + "ATOM 12531 HH12 ARG D 474 45.558 -39.551 -49.936 1.00 15.00 H \n" ++ + "ENDMDL \n" ++ + "trash strings 4\n" ++ + "MODEL 6 \n" ++ + "ATOM 12532 HH21 ARG D 474 47.457 -38.007 -47.445 1.00 15.00 H \n" ++ + "TER 12534 ARG D 474 \n" ++ + "ATOM 12533 HH22 ARG D 474 47.405 -39.268 -48.629 1.00 15.00 H \n" ++ + "HETATM12535 C1 NAG B 475 5.791 -20.194 -7.051 1.00 34.66 C \n" ++ + "HETATM12538 C4 NAG B 475 6.943 -19.507 -9.597 1.00 25.87 C \n" ++ + "CONECT conect conect\n" ++ + "ENDMDL\n" ++ + "MASTER 1 2 3 4 5 6 7 8\n" ++ + "trash strings 5\n" + ) + mt `shouldBe` Left "There are trash strings between model strings" + +pdbWithoutTrash :: PDB +pdbWithoutTrash = PDB { title = "STRUCTURE OF THE TRANSFORMED MONOCLINIC LYSOZYME" + , remarks = Data.Map.Strict.fromList [(Just 1, V.singleton "REFERENCE 1"), (Just 2, V.singleton "Reference 2_1/1")] + , models = V.fromList [V.fromList [V.fromList [Atom {atomSerial = 2033, atomName = " H ", atomAltLoc = ' ', atomResName = "CYS", + atomChainID = 'A', atomResSeq = 214, atomICode = ' ', atomX = -6.124, atomY = -27.225, + atomZ = 26.558, atomOccupancy = 1.0, atomTempFactor = 15.0, atomElement = " H", atomCharge = " "}] + ], + V.fromList [V.fromList [Atom {atomSerial = 2035, atomName = " N ", atomAltLoc = ' ', atomResName = "GLU", + atomChainID = 'B', atomResSeq = 1, atomICode = ' ', atomX = 18.637, atomY = -691.583, + atomZ = 66.852, atomOccupancy = 1.0, atomTempFactor = 118.48, atomElement = " N", atomCharge = " "}], + V.fromList [Atom {atomSerial = 12531, atomName = "HH12", atomAltLoc = ' ', atomResName = "ARG", + atomChainID = 'D', atomResSeq = 474, atomICode = ' ', atomX = 45.558, atomY = -39.551, atomZ = -49.936, + atomOccupancy = 1.0, atomTempFactor = 15.0, atomElement = " H", atomCharge = " "}] + ], + V.fromList [V.fromList [Atom {atomSerial = 12532, atomName = "HH21", atomAltLoc = ' ', atomResName = "ARG", + atomChainID = 'D', atomResSeq = 474, atomICode = ' ', atomX = 47.457, atomY = -38.007, atomZ = -47.445, + atomOccupancy = 1.0, atomTempFactor = 15.0, atomElement = " H", atomCharge = " "}, + Atom {atomSerial = 12533, atomName = "HH22", atomAltLoc = ' ', atomResName = "ARG", + atomChainID = 'D', atomResSeq = 474, atomICode = ' ', atomX = 47.405, atomY = -39.268, atomZ = -48.629, + atomOccupancy = 1.0, atomTempFactor = 15.0, atomElement = " H", atomCharge = " "}] + ] + ] + , otherFields = Data.Map.Strict.fromList [(HEADER, V.fromList [" header"]),(SEQRES,V.fromList [" seqres"]), + (CRYST1,V.fromList [" cryst1"]),(MASTER,V.fromList[" 1 2 3 4 5 6 7 8"])] + } + +onlyOneModelSpecP :: Spec +onlyOneModelSpecP = describe "Only One model." $ + it "correctly parses pdb with only one model without other field/title/trash/remarks strings" $ do + let mt = fromTextPDB . lenghtenLines $ T.pack ( "ATOM 2032 OXT CYS A 214 -4.546 -29.673 26.796 1.0 143.51 O \n" ++ + "ATOM 2033 H CYS A 214 -6.124 -27.225 26.558 1.00 15.00 H \n" ++ + "TER 2034 CYS A 214 \n" ++ + "ATOM 2035 N GLU B 1 18.637 -61.583 66.852 1.0 118.48 N \n" ++ + "TER 12534 ARG D 474 \n" ++ + "HETATM12535 C1 NAG B 475 5.791 -20.194 -7.051 1.00 34.66 C \n" ++ + "HETATM12538 C4 NAG B 475 6.943 -19.507 -9.597 1.00 25.87 C \n" ++ + "ATOM 2036 CA GLU B 1 19.722 -62.606 66.868 1.00 19.77 C \n" + ) + mt `shouldBe` Right ([], onlyOneModelPDB) + +onlyOneModelPDB :: PDB +onlyOneModelPDB = PDB { title = "" + , remarks = Data.Map.Strict.empty + , models = V.singleton $ V.fromList [ + V.fromList [ + Atom {atomSerial = 2032, atomName = " OXT", atomAltLoc = ' ', atomResName = "CYS", atomChainID = 'A', atomResSeq = 214, atomICode = ' ', + atomX = -4.546, atomY = -29.673, atomZ = 26.796, atomOccupancy = 1.0, atomTempFactor = 143.51, atomElement = " O", atomCharge = " "}, + Atom {atomSerial = 2033, atomName = " H ", atomAltLoc = ' ', atomResName = "CYS", atomChainID = 'A', atomResSeq = 214, atomICode = ' ', + atomX = -6.124, atomY = -27.225, atomZ = 26.558, atomOccupancy = 1.0, atomTempFactor = 15.0, atomElement = " H", atomCharge = " "}], + V.fromList [ + Atom {atomSerial = 2035, atomName = " N ", atomAltLoc = ' ', atomResName = "GLU", atomChainID = 'B', atomResSeq = 1, atomICode = ' ', + atomX = 18.637, atomY = -61.583, atomZ = 66.852, atomOccupancy = 1.0, atomTempFactor = 118.48, atomElement = " N", atomCharge = " "}, + Atom {atomSerial = 12535, atomName = " C1 ", atomAltLoc = ' ', atomResName = "NAG", atomChainID = 'B', atomResSeq = 475, atomICode = ' ', + atomX = 5.791, atomY = -20.194, atomZ = -7.051, atomOccupancy = 1.0, atomTempFactor = 34.66, atomElement = " C", atomCharge = " "}, + Atom {atomSerial = 12538, atomName = " C4 ", atomAltLoc = ' ', atomResName = "NAG", atomChainID = 'B', atomResSeq = 475, atomICode = ' ', + atomX = 6.943, atomY = -19.507, atomZ = -9.597, atomOccupancy = 1.0, atomTempFactor = 25.87, atomElement = " C", atomCharge = " "}, + Atom {atomSerial = 2036, atomName = " CA ", atomAltLoc = ' ', atomResName = "GLU", atomChainID = 'B', atomResSeq = 1, atomICode = ' ', + atomX = 19.722, atomY = -62.606, atomZ = 66.868, atomOccupancy = 1.0, atomTempFactor = 19.77, atomElement = " C", atomCharge = " "}] + ] + , otherFields = Data.Map.Strict.empty + } + +repeatedStringsSpecP :: Spec +repeatedStringsSpecP = describe "PDB with repeated other field strings." $ + it "correctly parses pdb with repeated other field strings (SOURCE)" $ do + let mt = fromTextPDB . lenghtenLines $ T.pack ( "HEADER header\n" ++ + "TITLE STRUCTURE OF THE TRANSFORMED MONOCLINIC LYSOZYME\n" ++ + "COMPND compnd\n" ++ + "EXPDTA expdta\n" ++ + "AUTHOR Masha\n" ++ + "SOURCE source1\n" ++ + "KEYWDS keywds\n" ++ + "REVDAT revdat\n" ++ + "REMARK 1 REFERENCE 1\n" ++ + "SEQRES seqres\n" ++ + "CRYST1 cryst1\n" ++ + "ORIGX1 origx1 n=1\n" ++ + "SOURCE source2\n" ++ + "SCALE2 sclaen n=2\n" ++ + "MASTER 1 2 3 4 5 6 7 8\n" ++ + "SOURCE source3\n" ++ + "END" + ) + mt `shouldBe` Right ([], repeatedStringsPDB) + +repeatedStringsPDB :: PDB +repeatedStringsPDB = PDB { title = "STRUCTURE OF THE TRANSFORMED MONOCLINIC LYSOZYME" + , remarks = Data.Map.Strict.singleton (Just 1) (V.singleton "REFERENCE 1") + , models = V.empty + , otherFields = Data.Map.Strict.fromList [(HEADER, V.fromList [" header"]),(COMPND,V.fromList[" compnd"]),(SOURCE,V.fromList[" source1", " source2", " source3"]),(KEYWDS,V.fromList[" keywds"]),(EXPDTA,V.fromList[" expdta"]),(AUTHOR,V.fromList[" Masha"]),(REVDAT,V.fromList[" revdat"]),(SEQRES,V.fromList [" seqres"]), + (CRYST1,V.fromList [" cryst1"]),(ORIGX1,V.fromList[" origx1 n=1"]),(SCALE2,V.fromList [" sclaen n=2"]),(MASTER,V.fromList[" 1 2 3 4 5 6 7 8"])] + } + +emptyRemarkSpecP :: Spec +emptyRemarkSpecP = describe "PDB with repeated remark strings without code." $ + it "correctly parses pdb with repeated remark strings without code" $ do + let mt = fromTextPDB . lenghtenLines $ T.pack ( "REMARK 111 remark111_1/2\n" ++ + "COMPND compnd\n" ++ + "SOURCE source1\n" ++ + "KEYWDS keywds\n" ++ + "REMARK 111 remark111_2/2\n" ++ + "REVDAT revdat\n" ++ + "REMARK 2 remark2_1/1\n" ++ + "SOURCE source2\n" ++ + "REMARK empty remark 1/2\n" ++ + "END \n" ++ + "REMARK empty remark 2/2\n" ++ + "SCALE2 sclaen n=2\n" + ) + mt `shouldBe` Right ([], pdbWithEmptyRemarks) + +pdbWithEmptyRemarks :: PDB +pdbWithEmptyRemarks = PDB { title = "" + , remarks = Data.Map.Strict.fromList [(Just 2, V.fromList ["remark2_1/1"]), (Just 111, V.fromList ["remark111_1/2", "remark111_2/2"]), (Nothing, V.fromList ["empty remark 1/2", "empty remark 2/2"])] + , models = V.empty + , otherFields = Data.Map.Strict.fromList [(COMPND,V.fromList[" compnd"]),(SOURCE,V.fromList[" source1", " source2"]), + (KEYWDS,V.fromList[" keywds"]), (REVDAT,V.fromList[" revdat"]), (SCALE2,V.fromList [" sclaen n=2"])] + } + +emptyModelSpecP :: Spec +emptyModelSpecP = describe "PDB with one empty model." $ + it "correctly parses pdb with one model without strings inside" $ do + let mt = fromTextPDB . lenghtenLines $ T.pack ( "trash strings 1\n" ++ + "HEADER header\n" ++ + "TITLE STRUCTURE OF THE TRANSFORMED MONOCLINIC LYSOZYME\n" ++ + "REMARK 1 REFERENCE 1\n" ++ + "trash strings 2\n" ++ + "REMARK 2 Reference 2_1/1\n" ++ + "MODEL 4 \n" ++ + "ENDMDL\n" ++ + "MASTER 1 2 3 4 5 6 7 8\n" ++ + "trash strings 5\n" + ) + mt `shouldBe` Right ([], pdbEmptyModel) + + +pdbEmptyModel :: PDB +pdbEmptyModel = PDB { title = "STRUCTURE OF THE TRANSFORMED MONOCLINIC LYSOZYME" + , remarks = Data.Map.Strict.fromList [(Just 1, V.singleton "REFERENCE 1"), (Just 2, V.singleton "Reference 2_1/1")] + , models = V.empty + , otherFields = Data.Map.Strict.fromList [(HEADER, V.fromList [" header"]), + (MASTER,V.fromList[" 1 2 3 4 5 6 7 8"])] + } + +lenghtenLines :: Text -> Text +lenghtenLines text = longLinedText + where + textLines = T.lines text + longTextLines = changeLine <$> textLines + desiredLength = 80 -- cause it is max length in standart pdb file + longLinedText = T.intercalate "\n" longTextLines + + changeLine :: Text -> Text + changeLine line | T.length line > desiredLength = T.take desiredLength line + | T.length line < desiredLength = line <> T.replicate spacesCount " " + | otherwise = line + where + spacesCount = desiredLength - T.length line + diff --git a/test/PDBSpec.hs b/test/PDBSpec.hs --- a/test/PDBSpec.hs +++ b/test/PDBSpec.hs @@ -1,428 +1,149 @@ {-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE RecordWildCards #-} module PDBSpec where -import Bio.PDB.Reader (fromTextPDB) -import Bio.PDB.Type (Atom (..), FieldType (..), PDB (..)) -import qualified Data.Map.Strict (empty, fromList, singleton) -import Data.Text as T (Text, intercalate, length, lines, pack, - replicate, take) -import qualified Data.Vector as V (empty, fromList, singleton) -import Test.Hspec +import Bio.PDB (modelsFromPDBFile) +import Bio.MAE (modelsFromMaeFile) +import Bio.Structure ( Model(..), Chain(..), Residue(..) + , Atom(..), Bond(..), GlobalID(..), LocalID(..)) +import Control.Monad.IO.Class (MonadIO, liftIO) +import Control.Exception (evaluate) +import Control.DeepSeq (force, NFData) -oneModelSpecP :: Spec -oneModelSpecP = describe "One model." $ - it "correctly parses pdb with only one model without strings \"MODEL\" & \"ENDMDL\"" $ do - let mt = fromTextPDB . lenghtenLines $ T.pack ( "HEADER header\n" ++ - "TITLE STRUCTURE OF THE TRANSFORMED MONOCLINIC LYSOZYME\n" ++ - "COMPND compnd\n" ++ - "SOURCE source\n" ++ - "KEYWDS keywds\n" ++ - "AUTHOR Masha\n" ++ - "REVDAT revdat\n" ++ - "REMARK 1 REFERENCE 1\n" ++ - "SEQRES seqres\n" ++ - "CRYST1 cryst1\n" ++ - "ORIGX1 origx1 n=1\n" ++ - "SCALE2 sclaen n=2\n" ++ - "ATOM 2032 OXT CYS A 214 -4.546 -29.673 26.796 1.0 143.51 O \n" ++ - "ATOM 2033 H CYS A 214 -6.124 -27.225 26.558 1.00 15.00 H \n" ++ - "TER 2034 CYS A 214 \n" ++ - "ATOM 2035 N GLU B 1 18.637 -61.583 66.852 1.0 118.48 N \n" ++ - "TER 12534 ARG D 474 \n" ++ - "HETATM12535 C1 NAG B 475 5.791 -20.194 -7.051 1.00 34.66 C \n" ++ - "HETATM12538 C4 NAG B 475 6.943 -19.507 -9.597 1.00 25.87 C \n" ++ - "ATOM 2036 CA GLU B 1 19.722 -62.606 66.868 1.00 19.77 C \n" ++ - "CONECT conect\n" ++ - "CONECT conect conect\n" ++ - "MASTER 1 2 3 4 5 6 7 8\n" ++ - "END" - ) - mt `shouldBe` Right ([], Right oneModelPDB) +import Data.Either (fromRight) +import Data.Vector (Vector) +import qualified Data.Vector as V (head, length, toList, concatMap) +import Data.List (find) +import Data.Text (Text) +import qualified Data.Text as T (pack) +import Data.Map.Strict (Map, (!)) +import qualified Data.Map.Strict as M (fromList) +import Data.Set (Set) +import qualified Data.Set as S (fromList, size, difference) -oneModelPDB :: PDB -oneModelPDB = PDB { title = "STRUCTURE OF THE TRANSFORMED MONOCLINIC LYSOZYME" - , remarks = Data.Map.Strict.singleton (Just 1) (V.singleton "REFERENCE 1") - , models = V.singleton $ V.fromList [V.fromList [Atom {atomSerial = 2032, atomName = " OXT", atomAltLoc = ' ', atomResName = "CYS", atomChainID = 'A', atomResSeq = 214, atomICode = ' ', - atomX = -4.546, atomY = -29.673, atomZ = 26.796, atomOccupancy = 1.0, atomTempFactor = 143.51, atomElement = " O", - atomCharge = " "}, - Atom {atomSerial = 2033, atomName = " H ", atomAltLoc = ' ', atomResName = "CYS", atomChainID = 'A', - atomResSeq = 214, atomICode = ' ', atomX = -6.124, atomY = -27.225, atomZ = 26.558, atomOccupancy = 1.0, atomTempFactor = 15.0, - atomElement = " H", atomCharge = " "}], - V.fromList [Atom {atomSerial = 2035, atomName = " N ", atomAltLoc = ' ', atomResName = "GLU", atomChainID = 'B', atomResSeq = 1, atomICode = ' ', - atomX = 18.637, atomY = -61.583, atomZ = 66.852, atomOccupancy = 1.0, atomTempFactor = 118.48, atomElement = " N", atomCharge = " "}, - Atom {atomSerial = 2036, atomName = " CA ", atomAltLoc = ' ', atomResName = "GLU", atomChainID = 'B', atomResSeq = 1, atomICode = ' ', - atomX = 19.722, atomY = -62.606, atomZ = 66.868, atomOccupancy = 1.0, atomTempFactor = 19.77, atomElement = " C", atomCharge = " "}] - ] - , otherFields = Data.Map.Strict.fromList [(HEADER, V.fromList [" header"]),(COMPND,V.fromList[" compnd"]),(SOURCE,V.fromList[" source"]),(KEYWDS,V.fromList[" keywds"]),(AUTHOR,V.fromList[" Masha"]),(REVDAT,V.fromList[" revdat"]),(SEQRES,V.fromList [" seqres"]), - (CRYST1,V.fromList [" cryst1"]),(ORIGX1,V.fromList[" origx1 n=1"]),(SCALE2,V.fromList [" sclaen n=2"]),(MASTER,V.fromList[" 1 2 3 4 5 6 7 8"])] - } +import Test.Hspec -manyModelsSpecP :: Spec -manyModelsSpecP = describe "Some models." $ - it "correctly parses pdb with many models - they have strings \"MODEL\" & \"ENDMDL\" and text has disordered strings" $ do - let mt = fromTextPDB . lenghtenLines $ T.pack ( "HEADER header\n" ++ - "TITLE STRUCTURE OF THE TRANSFORMED MONOCLINIC LYSOZYME\n" ++ - "COMPND compnd\n" ++ - "SOURCE source\n" ++ - "REMARK 1 REFERENCE 1\n" ++ - "REMARK 2 Reference 2_1/1\n" ++ - "CRYST1 cryst1\n" ++ - "ORIGX1 origx1 n=1\n" ++ - "SCALE2 sclaen n=2\n" ++ - "MODEL 4 \n" ++ - "ATOM 2033 H CYS A 214 -6.124 -27.225 26.558 1.00 15.00 H \n" ++ - "TER 2034 CYS A 214 \n" ++ - "ANISOU anisou\n" ++ - "ENDMDL \n" ++ - "MODEL 5 \n" ++ - "ATOM 2035 N GLU B 1 18.637-691.583 66.852 1.0 118.48 N \n" ++ - "ATOM 12531 HH12 ARG D 474 45.558 -39.551 -49.936 1.00 15.00 H \n" ++ - "ENDMDL \n" ++ - "MODEL 6 \n" ++ - "ATOM 12532 HH21 ARG D 474 47.457 -38.007 -47.445 1.00 15.00 H \n" ++ - "TER 12534 ARG D 474 \n" ++ - "ATOM 12533 HH22 ARG D 474 47.405 -39.268 -48.629 1.00 15.00 H \n" ++ - "HETATM12535 C1 NAG B 475 5.791 -20.194 -7.051 1.00 34.66 C \n" ++ - "HETATM12538 C4 NAG B 475 6.943 -19.507 -9.597 1.00 25.87 C \n" ++ - "CONECT conect conect\n" ++ - "CONECT conect\n" ++ - "ENDMDL\n" ++ - "SEQRES seqres\n" ++ - "KEYWDS keywds\n" ++ - "EXPDTA expdta\n" ++ - "AUTHOR Masha\n" ++ - "REVDAT revdat\n" ++ - "MASTER 1 2 3 4 5 6 7 8\n" - ) - mt `shouldBe` Right ([], Right manyModelsPDB) +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 -manyModelsPDB :: PDB -manyModelsPDB = PDB { title = "STRUCTURE OF THE TRANSFORMED MONOCLINIC LYSOZYME" - , remarks = Data.Map.Strict.fromList [(Just 1, V.singleton "REFERENCE 1"), (Just 2, V.singleton "Reference 2_1/1")] - , models = V.fromList [V.fromList [V.fromList [Atom {atomSerial = 2033, atomName = " H ", atomAltLoc = ' ', atomResName = "CYS", - atomChainID = 'A', atomResSeq = 214, atomICode = ' ', atomX = -6.124, atomY = -27.225, - atomZ = 26.558, atomOccupancy = 1.0, atomTempFactor = 15.0, atomElement = " H", atomCharge = " "}] - ], - V.fromList [V.fromList [Atom {atomSerial = 2035, atomName = " N ", atomAltLoc = ' ', atomResName = "GLU", - atomChainID = 'B', atomResSeq = 1, atomICode = ' ', atomX = 18.637, atomY = -691.583, - atomZ = 66.852, atomOccupancy = 1.0, atomTempFactor = 118.48, atomElement = " N", atomCharge = " "}], - V.fromList [Atom {atomSerial = 12531, atomName = "HH12", atomAltLoc = ' ', atomResName = "ARG", - atomChainID = 'D', atomResSeq = 474, atomICode = ' ', atomX = 45.558, atomY = -39.551, atomZ = -49.936, - atomOccupancy = 1.0, atomTempFactor = 15.0, atomElement = " H", atomCharge = " "}] - ], - V.fromList [V.fromList [Atom {atomSerial = 12532, atomName = "HH21", atomAltLoc = ' ', atomResName = "ARG", - atomChainID = 'D', atomResSeq = 474, atomICode = ' ', atomX = 47.457, atomY = -38.007, atomZ = -47.445, - atomOccupancy = 1.0, atomTempFactor = 15.0, atomElement = " H", atomCharge = " "}, - Atom {atomSerial = 12533, atomName = "HH22", atomAltLoc = ' ', atomResName = "ARG", - atomChainID = 'D', atomResSeq = 474, atomICode = ' ', atomX = 47.405, atomY = -39.268, atomZ = -48.629, - atomOccupancy = 1.0, atomTempFactor = 15.0, atomElement = " H", atomCharge = " "}] - ] - ] - , otherFields = Data.Map.Strict.fromList [(HEADER, V.fromList [" header"]),(COMPND,V.fromList[" compnd"]),(SOURCE,V.fromList[" source"]),(KEYWDS,V.fromList[" keywds"]),(EXPDTA,V.fromList[" expdta"]),(AUTHOR,V.fromList[" Masha"]),(REVDAT,V.fromList[" revdat"]),(SEQRES,V.fromList [" seqres"]), - (CRYST1,V.fromList [" cryst1"]),(ORIGX1,V.fromList[" origx1 n=1"]),(SCALE2,V.fromList [" sclaen n=2"]),(MASTER,V.fromList[" 1 2 3 4 5 6 7 8"])] - } + 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 -noModelsSpecP :: Spec -noModelsSpecP = describe "No models." $ - it "correctly parses pdb without models (no ATOM, TER, HETATM strings)" $ do - let mt = fromTextPDB . lenghtenLines $ T.pack ( "HEADER header\n" ++ - "TITLE STRUCTURE OF THE TRANSFORMED MONOCLINIC LYSOZYME\n" ++ - "COMPND compnd\n" ++ - "SOURCE source\n" ++ - "KEYWDS keywds\n" ++ - "EXPDTA expdta\n" ++ - "AUTHOR Masha\n" ++ - "REVDAT revdat\n" ++ - "REMARK 1 REFERENCE 1\n" ++ - "SEQRES seqres\n" ++ - "CRYST1 cryst1\n" ++ - "ORIGX1 origx1 n=1\n" ++ - "SCALE2 sclaen n=2\n" ++ - "MASTER 1 2 3 4 5 6 7 8\n" - ) - mt `shouldBe` Right ([], Right noModelsPDB) +-- 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", + "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" -noModelsPDB :: PDB -noModelsPDB = PDB { title = "STRUCTURE OF THE TRANSFORMED MONOCLINIC LYSOZYME" - , remarks = Data.Map.Strict.singleton (Just 1) (V.singleton "REFERENCE 1") - , models = V.empty - , otherFields = Data.Map.Strict.fromList [(HEADER, V.fromList [" header"]),(COMPND,V.fromList[" compnd"]),(SOURCE,V.fromList[" source"]),(KEYWDS,V.fromList[" keywds"]),(EXPDTA,V.fromList[" expdta"]),(AUTHOR,V.fromList[" Masha"]),(REVDAT,V.fromList[" revdat"]),(SEQRES,V.fromList [" seqres"]), - (CRYST1,V.fromList [" cryst1"]),(ORIGX1,V.fromList[" origx1 n=1"]),(SCALE2,V.fromList [" sclaen n=2"]),(MASTER,V.fromList[" 1 2 3 4 5 6 7 8"])] - } + let pdbBondCount = V.length $ modelBonds modelFromPDB + let maeBondCount = V.length $ modelBonds modelFromMae -allFieldsModelSpecP :: Spec -allFieldsModelSpecP = describe "PDB with all strings." $ - it "correctly parses pdb with all types of string" $ do - let mt = fromTextPDB . lenghtenLines $ T.pack ( "HEADER header\n" ++ - "OBSLTE obslte\n" ++ - "TITLE STRUCTURE OF THE TRANSFORMED MONOCLINIC LYSOZYME\n" ++ - "SPLIT split\n" ++ - "CAVEAT caveat\n" ++ - "COMPND compnd\n" ++ - "SOURCE source\n" ++ - "KEYWDS keywds\n" ++ - "EXPDTA expdta1\n" ++ - "EXPDTA expdta2\n" ++ - "NUMMDL nummdl\n" ++ - "MDLTYP mdltyp mdltyp\n" ++ - "AUTHOR Masha\n" ++ - "REVDAT revdat\n" ++ - "SPRSDE sprsde\n" ++ - "JRNL jrnl\n" ++ - "REMARK 1 REFERENCE 1\n" ++ - "DBREF dbref\n" ++ - "DBREF1 dbref1\n" ++ - "DBREF2 dbref2_1\n" ++ - "DBREF2 dbref2_2\n" ++ - "SEQADV seqadv\n" ++ - "SEQRES seqres\n" ++ - "MODRES modres\n" ++ - "HET het\n" ++ - "HETNAM hetnam\n" ++ - "HETSYN hetsyn\n" ++ - "FORMUL of love\n" ++ - "HELIX helix\n" ++ - "SHEET sheet\n" ++ - "SSBOND ssbond\n" ++ - "LINK link\n" ++ - "CISPEP cispep\n" ++ - "SITE site\n" ++ - "CRYST1 cryst1\n" ++ - "ORIGX1 origx1 n=1\n" ++ - "SCALE2 sclaen n=2\n" ++ - "MTRIX3 matrixn n=3\n" ++ - "MODEL 1\n" ++ - "ATOM 12532 HH21 ARG D 474 47.457 -38.007 -47.445 1.00 15.00 H \n" ++ - "TER 12534 ARG D 474 \n" ++ - "HETATM12535 C1 NAG B 475 5.791 -20.194 -7.051 1.00 34.66 C \n" ++ - "CONECT conect\n" ++ - "CONECT conect conect\n" ++ - "ENDMDL\n" ++ - "MASTER 1 2 3 4 5 6 7 8\n" ++ - "END" - ) - mt `shouldBe` Right ([], Right pdbWithAllFields) + it (tripeptideName ++ " equal bond count in Mae and PDB") $ pdbBondCount `shouldBe` maeBondCount -pdbWithAllFields :: PDB -pdbWithAllFields = PDB { title = "STRUCTURE OF THE TRANSFORMED MONOCLINIC LYSOZYME" - , models = V.fromList [V.fromList [V.singleton Atom {atomSerial = 12532, atomName = "HH21", atomAltLoc = ' ', atomResName = "ARG", - atomChainID = 'D', atomResSeq = 474, atomICode = ' ', atomX = 47.457, atomY = -38.007, atomZ = -47.445, - atomOccupancy = 1.0, atomTempFactor = 15.0, atomElement = " H", atomCharge = " "}]] - , remarks = Data.Map.Strict.singleton (Just 1) (V.singleton "REFERENCE 1") - , otherFields = Data.Map.Strict.fromList [(HEADER, V.fromList [" header"]),(OBSLTE, V.fromList [" obslte"]), (SPLIT, V.fromList ["split"]), (CAVEAT, V.fromList [" caveat"]), - (COMPND, V.fromList [" compnd"]),(SOURCE, V.fromList [" source"]),(KEYWDS, V.fromList [" keywds"]), - (EXPDTA, V.fromList [" expdta1"," expdta2"]),(NUMMDL, V.fromList [" nummdl"]),(MDLTYP, V.fromList [" mdltyp mdltyp"]), - (AUTHOR, V.fromList [" Masha"]),(REVDAT, V.fromList [" revdat"]),(SPRSDE, V.fromList [" sprsde"]), - (JRNL, V.fromList [" jrnl"]),(DBREF, V.fromList ["dbref"]), (DBREF1, V.fromList [" dbref1"]), (DBREF2, V.fromList [" dbref2_1"," dbref2_2"]), - (SEQADV, V.fromList [" seqadv"]),(SEQRES, V.fromList [" seqres"]),(MODRES, V.fromList [" modres"]), - (HET, V.fromList [" het"]), (HETNAM, V.fromList [" hetnam"]), (HETSYN, V.fromList [" hetsyn"]),(FORMUL, V.fromList [" of love"]),(HELIX, V.fromList ["helix"]), - (SHEET, V.fromList ["sheet"]),(SSBOND, V.fromList [" ssbond"]),(LINK, V.fromList [" link"]), - (CISPEP, V.fromList [" cispep"]),(SITE,V.fromList [" site"]),(CRYST1, V.fromList [" cryst1"]),(MTRIX3, V.fromList [" matrixn n=3"]), - (ORIGX1, V.fromList [" origx1 n=1"]),(SCALE2, V.fromList [" sclaen n=2"]),(MASTER, V.fromList [" 1 2 3 4 5 6 7 8"])] - } +bondsRestoringBiggerMoleculesSpec :: SpecWith () +bondsRestoringBiggerMoleculesSpec = describe "Bonds should be restored correctly in bigger molecules" $ do + checkBiggerMolecule "3mxw_ab_b" + checkBiggerMolecule "1vfb_ab_b" + 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 -emptySpecP :: Spec -emptySpecP = describe "empty PDB." $ - it "correctly parses empty pdb" $ do - let mt = fromTextPDB "" - mt `shouldBe` Right ([], Right emptyPdb) + it (moleculeName ++ " no dublicate bonds") $ length (doubleBonds modelFromPDB) `shouldBe` 0 -emptyPdb :: PDB -emptyPdb = PDB { title = "" - , models = V.empty - , remarks = Data.Map.Strict.empty - , otherFields = Data.Map.Strict.empty - } + let _globalBondSetPDB = globalBondSet modelFromPDB + let _globalBondSetMae = globalBondSet modelFromMae + let diffMaePDBGlobal = S.difference _globalBondSetMae _globalBondSetPDB + let diffPDBMaeGlobal = S.difference _globalBondSetPDB _globalBondSetMae -trashBetweenModelsSpecP :: Spec -trashBetweenModelsSpecP = describe "PDB has trash." $ - it "correctly parses pdb with trash string between models and other field strings" $ do - let mt = fromTextPDB . lenghtenLines $ T.pack ( "trash strings 1\n" ++ - "HEADER header\n" ++ - "TITLE STRUCTURE OF THE TRANSFORMED MONOCLINIC LYSOZYME\n" ++ - "REMARK 1 REFERENCE 1\n" ++ - "trash strings 2\n" ++ - "REMARK 2 Reference 2_1/1\n" ++ - "MODEL 4 \n" ++ - "ATOM 2033 H CYS A 214 -6.124 -27.225 26.558 1.00 15.00 H \n" ++ - "TER 2034 CYS A 214 \n" ++ - "ANISOU anisou\n" ++ - "ENDMDL \n" ++ - "CRYST1 cryst1\n" ++ - "trash strings 3\n" ++ - "SEQRES seqres\n" ++ - "MODEL 5 \n" ++ - "ATOM 2035 N GLU B 1 18.637-691.583 66.852 1.0 118.48 N \n" ++ - "ATOM 12531 HH12 ARG D 474 45.558 -39.551 -49.936 1.00 15.00 H \n" ++ - "ENDMDL \n" ++ - "trash strings 4\n" ++ - "MODEL 6 \n" ++ - "ATOM 12532 HH21 ARG D 474 47.457 -38.007 -47.445 1.00 15.00 H \n" ++ - "TER 12534 ARG D 474 \n" ++ - "ATOM 12533 HH22 ARG D 474 47.405 -39.268 -48.629 1.00 15.00 H \n" ++ - "HETATM12535 C1 NAG B 475 5.791 -20.194 -7.051 1.00 34.66 C \n" ++ - "HETATM12538 C4 NAG B 475 6.943 -19.507 -9.597 1.00 25.87 C \n" ++ - "CONECT conect conect\n" ++ - "ENDMDL\n" ++ - "MASTER 1 2 3 4 5 6 7 8\n" ++ - "trash strings 5\n" - ) - mt `shouldBe` Left "There are trash strings between model strings" + 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 -pdbWithoutTrash :: PDB -pdbWithoutTrash = PDB { title = "STRUCTURE OF THE TRANSFORMED MONOCLINIC LYSOZYME" - , remarks = Data.Map.Strict.fromList [(Just 1, V.singleton "REFERENCE 1"), (Just 2, V.singleton "Reference 2_1/1")] - , models = V.fromList [V.fromList [V.fromList [Atom {atomSerial = 2033, atomName = " H ", atomAltLoc = ' ', atomResName = "CYS", - atomChainID = 'A', atomResSeq = 214, atomICode = ' ', atomX = -6.124, atomY = -27.225, - atomZ = 26.558, atomOccupancy = 1.0, atomTempFactor = 15.0, atomElement = " H", atomCharge = " "}] - ], - V.fromList [V.fromList [Atom {atomSerial = 2035, atomName = " N ", atomAltLoc = ' ', atomResName = "GLU", - atomChainID = 'B', atomResSeq = 1, atomICode = ' ', atomX = 18.637, atomY = -691.583, - atomZ = 66.852, atomOccupancy = 1.0, atomTempFactor = 118.48, atomElement = " N", atomCharge = " "}], - V.fromList [Atom {atomSerial = 12531, atomName = "HH12", atomAltLoc = ' ', atomResName = "ARG", - atomChainID = 'D', atomResSeq = 474, atomICode = ' ', atomX = 45.558, atomY = -39.551, atomZ = -49.936, - atomOccupancy = 1.0, atomTempFactor = 15.0, atomElement = " H", atomCharge = " "}] - ], - V.fromList [V.fromList [Atom {atomSerial = 12532, atomName = "HH21", atomAltLoc = ' ', atomResName = "ARG", - atomChainID = 'D', atomResSeq = 474, atomICode = ' ', atomX = 47.457, atomY = -38.007, atomZ = -47.445, - atomOccupancy = 1.0, atomTempFactor = 15.0, atomElement = " H", atomCharge = " "}, - Atom {atomSerial = 12533, atomName = "HH22", atomAltLoc = ' ', atomResName = "ARG", - atomChainID = 'D', atomResSeq = 474, atomICode = ' ', atomX = 47.405, atomY = -39.268, atomZ = -48.629, - atomOccupancy = 1.0, atomTempFactor = 15.0, atomElement = " H", atomCharge = " "}] - ] - ] - , otherFields = Data.Map.Strict.fromList [(HEADER, V.fromList [" header"]),(SEQRES,V.fromList [" seqres"]), - (CRYST1,V.fromList [" cryst1"]),(MASTER,V.fromList[" 1 2 3 4 5 6 7 8"])] - } + let _localBondSetPDB = localBondSet modelFromPDB + let _localBondSetMae = localBondSet modelFromMae + let diffMaePDBLocal = S.difference _localBondSetMae _localBondSetPDB + let diffPDBMaeLocal = S.difference _localBondSetPDB _localBondSetMae -onlyOneModelSpecP :: Spec -onlyOneModelSpecP = describe "Only One model." $ - it "correctly parses pdb with only one model without other field/title/trash/remarks strings" $ do - let mt = fromTextPDB . lenghtenLines $ T.pack ( "ATOM 2032 OXT CYS A 214 -4.546 -29.673 26.796 1.0 143.51 O \n" ++ - "ATOM 2033 H CYS A 214 -6.124 -27.225 26.558 1.00 15.00 H \n" ++ - "TER 2034 CYS A 214 \n" ++ - "ATOM 2035 N GLU B 1 18.637 -61.583 66.852 1.0 118.48 N \n" ++ - "TER 12534 ARG D 474 \n" ++ - "HETATM12535 C1 NAG B 475 5.791 -20.194 -7.051 1.00 34.66 C \n" ++ - "HETATM12538 C4 NAG B 475 6.943 -19.507 -9.597 1.00 25.87 C \n" ++ - "ATOM 2036 CA GLU B 1 19.722 -62.606 66.868 1.00 19.77 C \n" - ) - mt `shouldBe` Right ([], Right onlyOneModelPDB) + 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 -onlyOneModelPDB :: PDB -onlyOneModelPDB = PDB { title = "" - , remarks = Data.Map.Strict.empty - , models = V.singleton $ V.fromList [V.fromList [Atom {atomSerial = 2032, atomName = " OXT", atomAltLoc = ' ', atomResName = "CYS", atomChainID = 'A', atomResSeq = 214, atomICode = ' ', - atomX = -4.546, atomY = -29.673, atomZ = 26.796, atomOccupancy = 1.0, atomTempFactor = 143.51, atomElement = " O", - atomCharge = " "}, - Atom {atomSerial = 2033, atomName = " H ", atomAltLoc = ' ', atomResName = "CYS", atomChainID = 'A', - atomResSeq = 214, atomICode = ' ', atomX = -6.124, atomY = -27.225, atomZ = 26.558, atomOccupancy = 1.0, atomTempFactor = 15.0, - atomElement = " H", atomCharge = " "}], - V.fromList [Atom {atomSerial = 2035, atomName = " N ", atomAltLoc = ' ', atomResName = "GLU", atomChainID = 'B', atomResSeq = 1, atomICode = ' ', - atomX = 18.637, atomY = -61.583, atomZ = 66.852, atomOccupancy = 1.0, atomTempFactor = 118.48, atomElement = " N", atomCharge = " "}, - Atom {atomSerial = 2036, atomName = " CA ", atomAltLoc = ' ', atomResName = "GLU", atomChainID = 'B', atomResSeq = 1, atomICode = ' ', - atomX = 19.722, atomY = -62.606, atomZ = 66.868, atomOccupancy = 1.0, atomTempFactor = 19.77, atomElement = " C", atomCharge = " "}] - ] - , otherFields = Data.Map.Strict.empty - } + localBondSet :: Model -> Set (Text, Int, Int, Int) -- (ChainID, ResidueNumber, LocalFrom, LocalTo) + localBondSet Model{..} = S.fromList $ do + Chain{..} <- V.toList modelChains + Residue{..} <- V.toList chainResidues + Bond (LocalID from) (LocalID to) _ <- V.toList resBonds -repeatedStringsSpecP :: Spec -repeatedStringsSpecP = describe "PDB with repeated other field strings." $ - it "correctly parses pdb with repeated other field strings (SOURCE)" $ do - let mt = fromTextPDB . lenghtenLines $ T.pack ( "HEADER header\n" ++ - "TITLE STRUCTURE OF THE TRANSFORMED MONOCLINIC LYSOZYME\n" ++ - "COMPND compnd\n" ++ - "EXPDTA expdta\n" ++ - "AUTHOR Masha\n" ++ - "SOURCE source1\n" ++ - "KEYWDS keywds\n" ++ - "REVDAT revdat\n" ++ - "REMARK 1 REFERENCE 1\n" ++ - "SEQRES seqres\n" ++ - "CRYST1 cryst1\n" ++ - "ORIGX1 origx1 n=1\n" ++ - "SOURCE source2\n" ++ - "SCALE2 sclaen n=2\n" ++ - "MASTER 1 2 3 4 5 6 7 8\n" ++ - "SOURCE source3\n" ++ - "END" - ) - mt `shouldBe` Right ([], Right repeatedStringsPDB) + [(chainName, resNumber, from, to), (chainName, resNumber, to, from)] -repeatedStringsPDB :: PDB -repeatedStringsPDB = PDB { title = "STRUCTURE OF THE TRANSFORMED MONOCLINIC LYSOZYME" - , remarks = Data.Map.Strict.singleton (Just 1) (V.singleton "REFERENCE 1") - , models = V.empty - , otherFields = Data.Map.Strict.fromList [(HEADER, V.fromList [" header"]),(COMPND,V.fromList[" compnd"]),(SOURCE,V.fromList[" source1", " source2", " source3"]),(KEYWDS,V.fromList[" keywds"]),(EXPDTA,V.fromList[" expdta"]),(AUTHOR,V.fromList[" Masha"]),(REVDAT,V.fromList[" revdat"]),(SEQRES,V.fromList [" seqres"]), - (CRYST1,V.fromList [" cryst1"]),(ORIGX1,V.fromList[" origx1 n=1"]),(SCALE2,V.fromList [" sclaen n=2"]),(MASTER,V.fromList[" 1 2 3 4 5 6 7 8"])] - } + doubleBonds :: Model -> [Bond GlobalID] + doubleBonds Model{..} = doubleBonds' (V.toList modelBonds) [] + where + doubleBonds' :: [Bond GlobalID] -> [Bond GlobalID] -> [Bond GlobalID] + doubleBonds' [] acc = acc + doubleBonds' (b:bs) acc = doubleBonds' bs . maybe acc (:acc) $ find (bondsEqual b) bs + bondsEqual :: Bond GlobalID -> Bond GlobalID -> Bool + bondsEqual b1 b2 = bondStart b1 == bondStart b2 && bondEnd b1 == bondEnd b2 -emptyRemarkSpecP :: Spec -emptyRemarkSpecP = describe "PDB with repeated remark strings without code." $ - it "correctly parses pdb with repeated remark strings without code" $ do - let mt = fromTextPDB . lenghtenLines $ T.pack ( "REMARK 111 remark111_1/2\n" ++ - "COMPND compnd\n" ++ - "SOURCE source1\n" ++ - "KEYWDS keywds\n" ++ - "REMARK 111 remark111_2/2\n" ++ - "REVDAT revdat\n" ++ - "REMARK 2 remark2_1/1\n" ++ - "SOURCE source2\n" ++ - "REMARK empty remark 1/2\n" ++ - "END \n" ++ - "REMARK empty remark 2/2\n" ++ - "SCALE2 sclaen n=2\n" - ) - mt `shouldBe` Right ([], Right pdbWithEmptyRemarks) + globalBondSet :: Model -> Set (Text,Text) + globalBondSet Model{..} = bondSet getGlobalID (chainsAtomMap modelChains) modelBonds -pdbWithEmptyRemarks :: PDB -pdbWithEmptyRemarks = PDB { title = "" - , remarks = Data.Map.Strict.fromList [(Just 2, V.fromList ["remark2_1/1"]), (Just 111, V.fromList ["remark111_1/2", "remark111_2/2"]), (Nothing, V.fromList ["empty remark 1/2", "empty remark 2/2"])] - , models = V.empty - , otherFields = Data.Map.Strict.fromList [(COMPND,V.fromList[" compnd"]),(SOURCE,V.fromList[" source1", " source2"]), - (KEYWDS,V.fromList[" keywds"]), (REVDAT,V.fromList[" revdat"]), (SCALE2,V.fromList [" sclaen n=2"])] - } + chainsAtomMap :: Vector Chain -> Map Int (Text, Atom) + chainsAtomMap chains = M.fromList $ concatMap chainAtomPreMap chains -emptyModelSpecP :: Spec -emptyModelSpecP = describe "PDB with one empty model." $ - it "correctly parses pdb with one model without strings inside" $ do - let mt = fromTextPDB . lenghtenLines $ T.pack ( "trash strings 1\n" ++ - "HEADER header\n" ++ - "TITLE STRUCTURE OF THE TRANSFORMED MONOCLINIC LYSOZYME\n" ++ - "REMARK 1 REFERENCE 1\n" ++ - "trash strings 2\n" ++ - "REMARK 2 Reference 2_1/1\n" ++ - "MODEL 4 \n" ++ - "ENDMDL\n" ++ - "MASTER 1 2 3 4 5 6 7 8\n" ++ - "trash strings 5\n" - ) - mt `shouldBe` Right ([], Right pdbEmptyModel) + chainAtomPreMap :: Chain -> [(Int, (Text, Atom))] + chainAtomPreMap Chain{..} = V.toList . fmap (\a -> (getGlobalID $ atomId a, (chainName, a))) $ V.concatMap resAtoms chainResidues + bondSet :: (a -> Int) -> Map Int (Text, Atom) -> Vector (Bond a) -> Set (Text,Text) + bondSet getID atomMap bonds = S.fromList $ do + Bond{..} <- V.toList bonds + let atomFromId = formAtomId $ atomMap ! getID bondStart + let atomToId = formAtomId $ atomMap ! getID bondEnd + [(atomFromId, atomToId), (atomToId, atomFromId)] + where + formAtomId :: (Text, Atom) -> Text + formAtomId (chainId, Atom{..}) = chainId <> "_" <> atomName <> "_" <> T.pack (show $ getGlobalID atomId) -pdbEmptyModel :: PDB -pdbEmptyModel = PDB { title = "STRUCTURE OF THE TRANSFORMED MONOCLINIC LYSOZYME" - , remarks = Data.Map.Strict.fromList [(Just 1, V.singleton "REFERENCE 1"), (Just 2, V.singleton "Reference 2_1/1")] - , models = V.empty - , otherFields = Data.Map.Strict.fromList [(HEADER, V.fromList [" header"]), - (MASTER,V.fromList[" 1 2 3 4 5 6 7 8"])] - } -lenghtenLines :: Text -> Text -lenghtenLines text = longLinedText - where - textLines = T.lines text - longTextLines = changeLine <$> textLines - desiredLength = 80 -- cause it is max length in standart pdb file - longLinedText = T.intercalate "\n" longTextLines +getStats :: Model -> (Int, Int, Int, Int) +getStats model = (globalBondCount, localBondCount, chainCount, atomCount) + where + globalBondCount = V.length $ modelBonds model + localBondCount = V.length . V.concatMap (V.concatMap resBonds . chainResidues) $ modelChains model + chainCount = V.length $ modelChains model + atomCount = V.length . V.concatMap (V.concatMap resAtoms . chainResidues) $ modelChains model - changeLine :: Text -> Text - changeLine line | T.length line > desiredLength = T.take desiredLength line - | T.length line < desiredLength = line <> T.replicate spacesCount " " - | otherwise = line - where - spacesCount = desiredLength - T.length line +ef :: NFData a => a -> IO a +ef = evaluate . force +firstPDBModel :: (MonadIO m) => FilePath -> m Model +firstPDBModel filepath = do + eitherPDB <- modelsFromPDBFile filepath + let (_, models) = fromRight undefined eitherPDB + liftIO . ef $ V.head models + +firstMaeModel :: (MonadIO m) => FilePath -> m Model +firstMaeModel filepath = do + eitherMae <- modelsFromMaeFile filepath + -- `evaluate . force` fails for some reason + pure . V.head $ fromRight undefined eitherMae diff --git a/test/Spec.hs b/test/Spec.hs --- a/test/Spec.hs +++ b/test/Spec.hs @@ -13,6 +13,7 @@ import System.IO import Test.Hspec import UniprotSpec +import PDBParserSpec main :: IO () main = do @@ -53,5 +54,8 @@ repeatedStringsSpecP emptyRemarkSpecP emptyModelSpecP + rawPDBToModelConversionSingleChainSpec + bondsRestoringTripeptideSpec + bondsRestoringBiggerMoleculesSpec -- Structure structureSpec