diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -2,6 +2,10 @@ ## [Unreleased] +## [0.1.3.18] - 2021-03-09 +### Fixed +- Lowercase insertion code parsing in PDB. + ## [0.1.3.17] - 2021-02-18 ### Fixed - PDB parsing. 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: 17226c38dac49205091506cf77e5ab4bff45f9ec44f1882a95c33e15cdb47e4e +-- hash: fe5f3bd32c503a44c34349ccfa5bf4ef5ec52540df5c2ce4637127a0b33394df name: cobot-io -version: 0.1.3.17 +version: 0.1.3.18 synopsis: Biological data file formats and IO description: Please see the README on GitHub at category: Bio @@ -69,6 +69,8 @@ Bio.Uniprot Bio.Uniprot.Parser Bio.Uniprot.Type + Bio.Utils.Map + Bio.Utils.Vector other-modules: Paths_cobot_io hs-source-dirs: diff --git a/src/Bio/MAE.hs b/src/Bio/MAE.hs --- a/src/Bio/MAE.hs +++ b/src/Bio/MAE.hs @@ -21,6 +21,7 @@ Model (..), Residue (..), SecondaryStructure (..), StructureModels (..)) +import qualified Bio.Utils.Map as M ((!?!)) import Control.Monad (join) import Control.Monad.IO.Class (MonadIO, liftIO) import Data.Attoparsec.Text (parseOnly) @@ -28,7 +29,7 @@ import Data.Function (on) import qualified Data.List as L (find, groupBy, sortOn) import Data.Map.Strict (Map) -import qualified Data.Map.Strict as M (fromList, lookup, (!)) +import qualified Data.Map.Strict as M (fromList, lookup) import Data.Maybe (catMaybes, fromJust) import Data.Text (Text) import qualified Data.Text as T (head, init, last, null, pack, @@ -62,7 +63,7 @@ modelsOf Mae{..} = V.fromList $ fmap blockToModel blocks where unsafeGetFromContentsMap :: FromMaeValue a => Map Text [MaeValue] -> Text -> Int -> a - unsafeGetFromContentsMap m name i = unsafeFromMaeValue $ (m M.! name) !! i + unsafeGetFromContentsMap m name i = unsafeFromMaeValue $ (m M.!?! name) !! i getFromContentsMap :: FromMaeValue a => Map Text [MaeValue] -> Text -> Int -> Maybe a getFromContentsMap m name i = join $ fromMaeValue . (!! i) <$> name `M.lookup` m @@ -71,7 +72,7 @@ blockToModel Block{..} = Model (atomsTableToChains atomsTable) bonds where atomsTable = findTable "m_atom" - numberOfAtoms = length $ atomsTable M.! "r_m_x_coord" + numberOfAtoms = length $ atomsTable M.!?! "r_m_x_coord" bondsTable = findTable "m_bond" (bondGraph, bonds) = bondsTableToGlobalBonds bondsTable @@ -89,7 +90,7 @@ bondsTableToGlobalBonds :: Map Text [MaeValue] -> (Map Int [(Int, Int)], Vector (Bond GlobalID)) bondsTableToGlobalBonds m = bimap toMap V.fromList bonds' where - numberOfBonds = length $ m M.! "i_m_from" + numberOfBonds = length $ m M.!?! "i_m_from" bonds' = unzip $ fmap indexToBond [0 .. numberOfBonds - 1] toMap :: [(Int, (Int, Int))] -> Map Int [(Int, Int)] @@ -150,7 +151,7 @@ toLocalBond :: Int -> (Int, Int) -> [Bond LocalID] toLocalBond x (y, o) | y `elem` group = pure $ Bond (LocalID x) - (LocalID $ globalToLocal M.! y) + (LocalID $ globalToLocal M.!?! y) o | otherwise = [] @@ -161,7 +162,7 @@ indexToAtom i = Atom (GlobalID i) (i + 1) (stripQuotes $ getFromContentsI "s_m_pdb_atom_name") - (elIndToElement M.! getFromContentsI "i_m_atomic_number") + (elIndToElement M.!?! getFromContentsI "i_m_atomic_number") coords (getFromContents 0 "i_m_formal_charge" i) (getFromContents 0 "r_m_pdb_tfactor" i) diff --git a/src/Bio/PDB.hs b/src/Bio/PDB.hs --- a/src/Bio/PDB.hs +++ b/src/Bio/PDB.hs @@ -11,21 +11,35 @@ import Bio.PDB.Reader (PDBWarnings, fromTextPDB) import qualified Bio.PDB.Type as PDB import Bio.PDB.Writer (pdbToFile, pdbToText) -import Bio.Structure +import Bio.Structure (Residue(..), + Bond, + Atom(..), + LocalID, + GlobalID(GlobalID, getGlobalID), + SecondaryStructure(Undefined), + Chain(..), + Model(Model, modelChains), + StructureSerializable(..), + StructureModels(..)) +import qualified Bio.Utils.Map as M ((!?!)) +import qualified Bio.Utils.Vector as V ((!?!)) import Control.Arrow ((&&&)) import Control.Lens ((^.)) import Control.Monad (join) import Control.Monad.IO.Class (MonadIO, liftIO) import Data.List (sort) -import Data.Map (Map) -import qualified Data.Map as M (fromList, (!)) +import Data.Map.Strict (Map) +import qualified Data.Map.Strict as M (fromList) import Data.Maybe (fromMaybe) import Data.Text (Text) import qualified Data.Text as T (head, pack, singleton, strip, unpack) import Data.Text.IO as TIO (readFile) import Data.Vector (Vector) -import qualified Data.Vector as V +import qualified Data.Vector as V (toList, + concat, + fromList, + length) import Linear.V3 (V3 (..), _x, _y, _z) import Text.Read (readMaybe) @@ -33,9 +47,9 @@ modelsOf PDB.PDB {..} = fmap mkModel models where mkModel :: PDB.Model -> Model - mkModel model = case length atomToNilBasedIndex == length allModelAtoms of - False -> error "Mapping from PDB id to nil based index must be a bijection." - True -> Model (fmap mkChain model) (restoreModelGlobalBonds atomToNilBasedIndex model) + mkModel model = if length atomToNilBasedIndex == length allModelAtoms + then Model (fmap mkChain model) (restoreModelGlobalBonds atomToNilBasedIndex model) + else error "Mapping from PDB id to nil based index must be a bijection." where atomToNilBasedIndex :: Map PDB.Atom Int atomToNilBasedIndex = M.fromList $ allModelAtoms `zip` [0..] @@ -49,27 +63,27 @@ mkChainName :: PDB.Chain -> Text mkChainName = T.singleton . PDB.atomChainID . safeFirstAtom - mkChainResidues :: PDB.Chain -> V.Vector Residue + mkChainResidues :: PDB.Chain -> 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 + safeFirstAtom :: 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 :: Map Text (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) + (localBondsMap M.!?! residueID firstResidueAtom) Undefined -- now we do not read secondary structure "" -- chemical component type?! where firstResidueAtom = head atoms' mkAtom :: PDB.Atom -> Atom - mkAtom atom@PDB.Atom{..} = Atom (GlobalID $ atomToNilBasedIndex M.! atom) + mkAtom atom@PDB.Atom{..} = Atom (GlobalID $ atomToNilBasedIndex M.!?! atom) atomSerial (T.strip atomName) atomElement @@ -78,10 +92,10 @@ atomTempFactor atomOccupancy -modelsFromPDBFile :: (MonadIO m) => FilePath -> m (Either Text ([PDBWarnings], V.Vector Model)) +modelsFromPDBFile :: (MonadIO m) => FilePath -> m (Either Text ([PDBWarnings], Vector Model)) modelsFromPDBFile = liftIO . fmap modelsFromPDBText . TIO.readFile -modelsFromPDBText :: Text -> Either Text ([PDBWarnings], V.Vector Model) +modelsFromPDBText :: Text -> Either Text ([PDBWarnings], Vector Model) modelsFromPDBText pdbText = do (warnings, parsedPDB) <- fromTextPDB pdbText let models = modelsOf parsedPDB diff --git a/src/Bio/PDB/Functions.hs b/src/Bio/PDB/Functions.hs --- a/src/Bio/PDB/Functions.hs +++ b/src/Bio/PDB/Functions.hs @@ -2,19 +2,25 @@ ( 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) +import qualified Bio.PDB.Type as PDB (Atom (..)) +import qualified Bio.Utils.Map as M ((!?!)) +import Data.Map.Strict (Map) +import qualified Data.Map.Strict as M (fromList) +import Data.List (groupBy, + sortOn) +import Data.Vector (Vector) +import qualified Data.Vector as V (toList) +import Data.Char (toUpper) + 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) + sortOnResidue PDB.Atom{..} = atomSerial * 100 + (insertionCodeSortingCorrections M.!?! toUpper atomICode) + insertionCodeSortingCorrections :: Map Char Int insertionCodeSortingCorrections = M.fromList $ zip (' ':['A'..'Z']) [0..] diff --git a/src/Bio/Structure/Functions.hs b/src/Bio/Structure/Functions.hs --- a/src/Bio/Structure/Functions.hs +++ b/src/Bio/Structure/Functions.hs @@ -12,7 +12,8 @@ localBonds, residues) import Control.Lens (Traversal', each, (%~), (&)) import Data.Map.Strict (Map) -import qualified Data.Map.Strict as M (fromList, (!), (!?)) +import qualified Data.Map.Strict as M (fromList, (!?)) +import qualified Bio.Utils.Map as M ((!?!)) import Data.Set (Set) import qualified Data.Set as S (fromList, notMember, unions) import Data.Text (Text) @@ -84,8 +85,8 @@ leaveBond (Bond (LocalID l) (LocalID r) _) = l `notElem` indsToDelete && r `notElem` indsToDelete modifyBond :: Bond LocalID -> Bond LocalID - modifyBond (Bond (LocalID l) (LocalID r) t) = Bond (LocalID $ oldIndsToNew M.! l) - (LocalID $ oldIndsToNew M.! r) + modifyBond (Bond (LocalID l) (LocalID r) t) = Bond (LocalID $ oldIndsToNew M.!?! l) + (LocalID $ oldIndsToNew M.!?! r) t newInd :: Int -> Int diff --git a/src/Bio/Utils/Map.hs b/src/Bio/Utils/Map.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Utils/Map.hs @@ -0,0 +1,13 @@ +module Bio.Utils.Map + ( (!?!) + ) where + + +import GHC.Stack (HasCallStack) +import Data.Map.Strict (Map, + (!?)) +import Data.Maybe (fromMaybe) + +infix 9 !?! +(!?!) :: (HasCallStack, Ord k, Show k, Show a) => Map k a -> k -> a +(!?!) m k = fromMaybe (error $ "cobot-io: No key " ++ show k ++ " in Map: " ++ show m) $ m !? k diff --git a/src/Bio/Utils/Vector.hs b/src/Bio/Utils/Vector.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Utils/Vector.hs @@ -0,0 +1,17 @@ +module Bio.Utils.Vector + ( (!?!) + ) where + + +import GHC.Stack (HasCallStack) +import Data.Vector (Vector) +import qualified Data.Vector as V ((!?), + length) +import Data.Maybe (fromMaybe) + +infix 9 !?! +(!?!) :: (HasCallStack, Show a) => Vector a -> Int -> a +(!?!) v i = fromMaybe (error msg) $ v V.!? i + where + msg :: String + msg = "cobot-io: index " ++ show i ++ " is out of bounds. Vector length is : " ++ show (V.length v) diff --git a/test/PDBSpec.hs b/test/PDBSpec.hs --- a/test/PDBSpec.hs +++ b/test/PDBSpec.hs @@ -130,6 +130,12 @@ formAtomId :: (Text, Atom) -> Text formAtomId (chainId, Atom{..}) = chainId <> "_" <> atomName <> "_" <> T.pack (show $ getGlobalID atomId) +lowercaseInsertionCodeSpec :: SpecWith () +lowercaseInsertionCodeSpec = describe "Lowercase insertion code in PDB" $ do + it "should be processed normally" $ do + model <- firstPDBModel "test/PDB/6a0z.HLHL.pdb" >>= ef + + V.length (modelChains model) `shouldBe` 2 getStats :: Model -> (Int, Int, Int, Int) getStats model = (globalBondCount, localBondCount, chainCount, atomCount) diff --git a/test/Spec.hs b/test/Spec.hs --- a/test/Spec.hs +++ b/test/Spec.hs @@ -59,6 +59,7 @@ bondsRestoringTripeptideSpec bondsRestoringBiggerMoleculesSpec withSegmentIdentifierModelSpecP + lowercaseInsertionCodeSpec pdbWriterSpec -- Structure structureSpec