diff --git a/PDBtools.cabal b/PDBtools.cabal --- a/PDBtools.cabal +++ b/PDBtools.cabal @@ -1,5 +1,5 @@ Name: PDBtools -Version: 0.0.1 +Version: 0.0.2 License: GPL-3 License-file: LICENSE Cabal-Version: >= 1.6 @@ -31,8 +31,10 @@ containers exposed-modules: - PDBtools.PDBparse - PDBtools.PDButil + PDBtools.Base + PDBtools.Residues + PDButil.PDBparse + PDButil.Vectors ghc-options: diff --git a/PDBtools/Base.hs b/PDBtools/Base.hs new file mode 100644 --- /dev/null +++ b/PDBtools/Base.hs @@ -0,0 +1,149 @@ +-- Module : PDBtools +-- Copyright : (c) 2012 Grant Rotskoff +-- License : GPL-3 +-- +-- Maintainer : gmr1887@gmail.com +-- Stability : experimental + + +-- The suite of source files in the PDButil directory are meant to complement high-throughput analyses +-- of three-dimensional protein structure data in the PDB format. Because the source files rely heavily +-- on one another, it is convenient to import them all into a single module. Examples of analysis projects +-- are available at http://www.github.com/rotskoff/Haskell-PDB-Utilities + + +module PDBtools.Base where + +-- Long list of imports... +import PDButil.PDBparse +import PDButil.Vectors +import Data.List +import Data.ByteString.Char8 (ByteString) +import qualified Data.ByteString.Char8 as B +import qualified Data.Map as Map +import Data.Maybe + +-- Pull all atoms of a given name from a list of atoms +atomtype :: String -> [Atom] -> [Atom] +atomtype t atmlist = filter matcht atmlist where + matcht a = atype a == B.pack t + +-- Match the atom's name in the PDB file rather than the underlying type +atomname :: String -> [Atom] -> [Atom] +atomname t atmlist = filter matcht atmlist where + matcht a = name a == B.pack t + +-- Match the residue type, input the three letter abbreviation +restype :: String -> [Atom] -> [Atom] +restype t atmlist = filter matcht atmlist where + matcht a = resname a == B.pack t + +-- Extract the list of alpha-Carbons from a protein +backbone :: Protein -> [Atom] +backbone = atomname "CA" . atoms + +-- Extract the list of residue name, residue number pairs +resSeq :: Protein -> [(ByteString,Int)] +resSeq p = zip (map resname bbatms) (map resid bbatms) where + bbatms = backbone p + +--Naive homology calculation +--this could be greatly improved by some alignment effort +homology :: Protein -> Protein -> Double +homology a b = (fromIntegral identities) / (fromIntegral totalLength) where + identities = 2 * length (resSeq a `intersect` resSeq b) + totalLength = length(resSeq a) + length(resSeq b) + +-- Euclidean Distance between two atoms +distance :: Atom -> Atom -> Double +distance a1 a2 = norm $ vSub (coords a1) (coords a2) + +-- Root Mean Squared Deviation, a measure of the total distance change +-- Only well-defined if you input the same protein! +rmsd :: [Atom] -> [Atom] -> Double +rmsd atms atms' = sqrt $ avg sqdist where + avg ds = (1/fromIntegral(length(ds))) * sum(ds) + sqdist = map (^2) $ zipWith distance (atms) (atms') + +-- Collect all the atoms within a given distance +within :: Double -> Atom -> [Atom] -> [Atom] +within range a = (delete a) . filter withinRange where + withinRange a' = (distance a a') <= range + +withinClusive :: Double -> Atom -> [Atom] -> [Atom] +withinClusive range a = filter withinRange where + withinRange a' = (distance a a') <= range + +-- Centers the list of atoms around the specified atom +center :: Atom -> [Atom] -> [Atom] +center a = a `shift` [0,0,0] + +-- Shift the entire atomlist by specifying the new location of a single atom +shift :: Atom -> [Double] -> [Atom] -> [Atom] +shift a newCoords as = map (\s -> s {coords = (translate s)}) as where + translate s = vAdd (coords s) shiftFactor + shiftFactor = vSub newCoords (coords a) + +-- Global translate by a vector +translateBy :: [Double] -> [Atom] -> [Atom] +translateBy vect = map (\s -> s {coords = (vAdd (coords s) vect)}) + +-- Compute the angle between three atoms, return value in radians! +atmAngle :: Atom -> Atom -> Atom -> Double +atmAngle a b c = angle baVec bcVec where + baVec = (coords a) `vSub` (coords b) + bcVec = (coords c) `vSub` (coords c) + +-- Convert a protein to FASTA sequence format +-- TODO, headers in FASTA file spec +protein2fasta :: Protein -> ByteString +protein2fasta protein = B.pack $ concatMap (\s -> convert (resname s)) (backbone protein) where + convert name = fromJust $ Map.lookup (B.unpack name) resMap + resMap = Map.fromList + [("ALA","A"), + ("CYS","C"), + ("ASP","D"), + ("GLU","E"), + ("PHE","F"), + ("GLY","G"), + ("HIS","H"), + ("ILE","I"), + ("LYS","K"), + ("LEU","L"), + ("MET","M"), + ("ASN","N"), + ("PYL","O"), + ("PRO","P"), + ("GLN","Q"), + ("ARG","R"), + ("SER","S"), + ("THR","T"), + --Selenocysteine + ("VAL","V"), + ("TRP","W"), + ("TYR","Y")] + -- otherwise, use 'X' + +{- +-- TODO Fix so that it works! +rotateAboutOrigin :: [Atom] -> Atom -> [Double] -> [Atom] +rotateAboutOrigin atms tracer destination = map (\s -> s {coords = (translate (coords s))}) atms where + [x,y,z] = destination + [x',y',z'] = coords tracer + psi = angle [y',z'] [y,z] --yz angle + p = angle [x',z'] [x,z] --xz angle + phi = angle [x',y'] [x,y] --xy angle + translate = vRotate3d theta psi phi + +-- Rotate a list of atoms about a fixed atom by moving the selected atom to a specified destination +rotate :: [Atom] -> Atom -> Atom -> [Double] -> [Atom] +rotate atms pivot tracer destination = translateBy (coords pivot) rotatedAtms where + centeredAtPivot = center pivot atms + rotatedAtms = rotateAboutOrigin centeredAtPivot tracer destination + +-} + + + + + diff --git a/PDBtools/PDBparse.hs b/PDBtools/PDBparse.hs deleted file mode 100644 --- a/PDBtools/PDBparse.hs +++ /dev/null @@ -1,83 +0,0 @@ -module PDBtools.PDBparse where - -import Data.ByteString.Char8 (ByteString) -import qualified Data.ByteString.Char8 as B -import System.IO (FilePath) - -data Atom = Atom { name :: ByteString, - atid :: Int, - chain :: ByteString, - resid :: Int, - resname :: ByteString, - coords :: [Double], - aField :: Double, - bField :: Double, - atype :: ByteString } - deriving (Show,Eq) - -data Protein = Protein { atoms :: [Atom] } - deriving (Show) - ---Sample record: --- ATOM 1 N ASP A 28 52.958 39.871 41.308 1.00 89.38 N - -{- We only want record lines that begin with ATOM and HETATM - ATOM lines contain the coordinates of the protein(s) in a PDB file - HETATM lines (short for heteroatom) contain coordinate information for - other molecules present in the structure... ligands, DNA, RNA, waters, etc. -} - -parseAtom :: ByteString -> Atom -parseAtom record = Atom { name = pull 13 16, - atid = rpull 7 11, - chain = pull 22 22, - resid = rpull 23 26, - resname = pull 18 20, - coords = [rpull 31 38,rpull 39 46,rpull 47 54], - aField = rpull 55 60, - bField = rpull 61 66, - atype = pull 77 78 } where - - --Hard coded parsing of the PDB record for coordinate types - --I've encountered this "repacking for comparison in expert code, - --but it seems like comparison should be possible some other way - - pull m n = cutspace $ B.drop (m-1) $ B.take n record - rpull m n = read $ B.unpack $ pull m n - cutspace = B.pack . filter (/=' ') . B.unpack - - -isAtom :: ByteString -> Bool -isAtom line = (B.take 4 line) == (B.pack "ATOM") - -isHETATM :: ByteString -> Bool -isHETATM line = (B.take 6 line) == (B.pack "HETATM") - - -parse :: FilePath -> IO ([Protein],[Atom]) -parse pdb = do - let input = B.readFile pdb - bstring <- input - let atms = map parseAtom $ filter isAtom (B.lines bstring) - let hetatms = map parseAtom $ filter isHETATM (B.lines bstring) - return (splitChains atms, hetatms) - -parseCofactorOnly :: FilePath -> IO [Atom] -parseCofactorOnly pdb = do - bstring <- B.readFile pdb - let hetatms = map parseAtom $ filter isHETATM (B.lines bstring) - return hetatms - -parseProteinOnly :: FilePath -> IO [Protein] -parseProteinOnly pdb = do - bstring <- B.readFile pdb - let atms = map parseAtom $ filter isAtom (B.lines bstring) - return $ splitChains atms - -splitChains :: [Atom] -> [Protein] -splitChains [] = [] -splitChains contents = [Protein {atoms = chain1}] ++ splitChains remainder where - chain1 = takeWhile (\s -> id == chain s) contents - remainder = dropWhile (\s -> id == chain s) contents - id = chain (head contents) - - diff --git a/PDBtools/PDButil.hs b/PDBtools/PDButil.hs deleted file mode 100644 --- a/PDBtools/PDButil.hs +++ /dev/null @@ -1,136 +0,0 @@ ---PDB Utilities, functions for basic manipulation of PDB data -module PDBtools.PDButil where - --- Long list of imports... -import PDBtools.PDBparse -import PDBtools.Vectors -import Data.List -import Data.ByteString.Char8 (ByteString) -import qualified Data.ByteString.Char8 as B -import qualified Data.Map as Map -import Data.Maybe - --- Pull all atoms of a given name from a list of atoms -atomtype :: String -> [Atom] -> [Atom] -atomtype t atmlist = filter matcht atmlist where - matcht a = atype a == B.pack t - --- Match the atom's name in the PDB file rather than the underlying type -atomname :: String -> [Atom] -> [Atom] -atomname t atmlist = filter matcht atmlist where - matcht a = name a == B.pack t - --- Match the residue type, input the three letter abbreviation -restype :: String -> [Atom] -> [Atom] -restype t atmlist = filter matcht atmlist where - matcht a = resname a == B.pack t - --- Extract the list of alpha-Carbons from a protein -backbone :: Protein -> [Atom] -backbone = atomname "CA" . atoms - --- Extract the list of residue name, residue number pairs -resSeq :: Protein -> [(ByteString,Int)] -resSeq p = zip (map resname bbatms) (map resid bbatms) where - bbatms = backbone p - ---Naive homology calculation ---this could be greatly improved by some alignment effort -homology :: Protein -> Protein -> Double -homology a b = (fromIntegral identities) / (fromIntegral totalLength) where - identities = 2 * length (resSeq a `intersect` resSeq b) - totalLength = length(resSeq a) + length(resSeq b) - --- Euclidean Distance between two atoms -distance :: Atom -> Atom -> Double -distance a1 a2 = norm $ vSub (coords a1) (coords a2) - --- Root Mean Squared Deviation, a measure of the total distance change --- Only well-defined if you input the same protein! -rmsd :: [Atom] -> [Atom] -> Double -rmsd atms atms' = sqrt $ avg sqdist where - avg ds = (1/fromIntegral(length(ds))) * sum(ds) - sqdist = map (^2) $ zipWith distance (atms) (atms') - --- Collect all the atoms within a given distance -within :: Double -> Atom -> [Atom] -> [Atom] -within range a = (delete a) . filter withinRange where - withinRange a' = (distance a a') <= range - -withinClusive :: Double -> Atom -> [Atom] -> [Atom] -withinClusive range a = filter withinRange where - withinRange a' = (distance a a') <= range - --- Centers the list of atoms around the specified atom -center :: Atom -> [Atom] -> [Atom] -center a = a `shift` [0,0,0] - --- Shift the entire atomlist by specifying the new location of a single atom -shift :: Atom -> [Double] -> [Atom] -> [Atom] -shift a newCoords as = map (\s -> s {coords = (translate s)}) as where - translate s = vAdd (coords s) shiftFactor - shiftFactor = vSub newCoords (coords a) - --- Global translate by a vector -translateBy :: [Double] -> [Atom] -> [Atom] -translateBy vect = map (\s -> s {coords = (vAdd (coords s) vect)}) - --- Compute the angle between three atoms, return value in radians! -atmAngle :: Atom -> Atom -> Atom -> Double -atmAngle a b c = angle baVec bcVec where - baVec = (coords a) `vSub` (coords b) - bcVec = (coords c) `vSub` (coords c) - --- Convert a protein to FASTA sequence format --- TODO, headers in FASTA file spec -protein2fasta :: Protein -> ByteString -protein2fasta protein = B.pack $ concatMap (\s -> convert (resname s)) (backbone protein) where - convert name = fromJust $ Map.lookup (B.unpack name) resMap - resMap = Map.fromList - [("ALA","A"), - ("CYS","C"), - ("ASP","D"), - ("GLU","E"), - ("PHE","F"), - ("GLY","G"), - ("HIS","H"), - ("ILE","I"), - ("LYS","K"), - ("LEU","L"), - ("MET","M"), - ("ASN","N"), - ("PYL","O"), - ("PRO","P"), - ("GLN","Q"), - ("ARG","R"), - ("SER","S"), - ("THR","T"), - --Selenocysteine - ("VAL","V"), - ("TRP","W"), - ("TYR","Y")] - -- otherwise, use 'X' - -{- --- TODO Fix so that it works! -rotateAboutOrigin :: [Atom] -> Atom -> [Double] -> [Atom] -rotateAboutOrigin atms tracer destination = map (\s -> s {coords = (translate (coords s))}) atms where - [x,y,z] = destination - [x',y',z'] = coords tracer - psi = angle [y',z'] [y,z] --yz angle - p = angle [x',z'] [x,z] --xz angle - phi = angle [x',y'] [x,y] --xy angle - translate = vRotate3d theta psi phi - --- Rotate a list of atoms about a fixed atom by moving the selected atom to a specified destination -rotate :: [Atom] -> Atom -> Atom -> [Double] -> [Atom] -rotate atms pivot tracer destination = translateBy (coords pivot) rotatedAtms where - centeredAtPivot = center pivot atms - rotatedAtms = rotateAboutOrigin centeredAtPivot tracer destination - --} - - - - - diff --git a/PDBtools/Residues.hs b/PDBtools/Residues.hs new file mode 100644 --- /dev/null +++ b/PDBtools/Residues.hs @@ -0,0 +1,34 @@ +-- Residues centered at the Carbon-Alpha, some sort of directionality constraint + +module PDBtools.Residues where + +import PDBtools.Base +import PDButil.PDBparse + +--One should really only use these methods on proteins, but for the sake of composing selections, the input form is [Atom] + +charged :: [Atom] -> [Atom] +charged atms = restype "ASP" atms + ++ restype "GLU" atms + ++ restype "ARG" atms + ++ restype "LYS" atms + ++ restype "HIS" atms + +uncharged :: [Atom] -> [Atom] +uncharged atms = filter (\s -> elem s $ charged atms) atms + +polar :: [Atom] -> [Atom] +polar atms = charged atms + ++ restype "SER" atms + ++ restype "THR" atms + ++ restype "ASN" atms + ++ restype "GLN" atms + +nonpolar :: [Atom] -> [Atom] +nonpolar atms = filter (\s -> elem s $ polar atms) atms + +hydrophobic :: [Atom] -> [Atom] +hydrophobic = nonpolar + +hydrophillic :: [Atom] -> [Atom] +hydrophillic = polar diff --git a/PDButil/PDBparse.hs b/PDButil/PDBparse.hs new file mode 100644 --- /dev/null +++ b/PDButil/PDBparse.hs @@ -0,0 +1,83 @@ +module PDButil.PDBparse where + +import Data.ByteString.Char8 (ByteString) +import qualified Data.ByteString.Char8 as B +import System.IO (FilePath) + +data Atom = Atom { name :: ByteString, + atid :: Int, + chain :: ByteString, + resid :: Int, + resname :: ByteString, + coords :: [Double], + aField :: Double, + bField :: Double, + atype :: ByteString } + deriving (Show,Eq) + +data Protein = Protein { atoms :: [Atom] } + deriving (Show) + +--Sample record: +-- ATOM 1 N ASP A 28 52.958 39.871 41.308 1.00 89.38 N + +{- We only want record lines that begin with ATOM and HETATM + ATOM lines contain the coordinates of the protein(s) in a PDB file + HETATM lines (short for heteroatom) contain coordinate information for + other molecules present in the structure... ligands, DNA, RNA, waters, etc. -} + +parseAtom :: ByteString -> Atom +parseAtom record = Atom { name = pull 13 16, + atid = rpull 7 11, + chain = pull 22 22, + resid = rpull 23 26, + resname = pull 18 20, + coords = [rpull 31 38,rpull 39 46,rpull 47 54], + aField = rpull 55 60, + bField = rpull 61 66, + atype = pull 77 78 } where + + --Hard coded parsing of the PDB record for coordinate types + --I've encountered this "repacking for comparison in expert code, + --but it seems like comparison should be possible some other way + + pull m n = cutspace $ B.drop (m-1) $ B.take n record + rpull m n = read $ B.unpack $ pull m n + cutspace = B.pack . filter (/=' ') . B.unpack + + +isAtom :: ByteString -> Bool +isAtom line = (B.take 4 line) == (B.pack "ATOM") + +isHETATM :: ByteString -> Bool +isHETATM line = (B.take 6 line) == (B.pack "HETATM") + + +parse :: FilePath -> IO ([Protein],[Atom]) +parse pdb = do + let input = B.readFile pdb + bstring <- input + let atms = map parseAtom $ filter isAtom (B.lines bstring) + let hetatms = map parseAtom $ filter isHETATM (B.lines bstring) + return (splitChains atms, hetatms) + +parseCofactorOnly :: FilePath -> IO [Atom] +parseCofactorOnly pdb = do + bstring <- B.readFile pdb + let hetatms = map parseAtom $ filter isHETATM (B.lines bstring) + return hetatms + +parseProteinOnly :: FilePath -> IO [Protein] +parseProteinOnly pdb = do + bstring <- B.readFile pdb + let atms = map parseAtom $ filter isAtom (B.lines bstring) + return $ splitChains atms + +splitChains :: [Atom] -> [Protein] +splitChains [] = [] +splitChains contents = [Protein {atoms = chain1}] ++ splitChains remainder where + chain1 = takeWhile (\s -> id == chain s) contents + remainder = dropWhile (\s -> id == chain s) contents + id = chain (head contents) + + diff --git a/PDButil/Vectors.hs b/PDButil/Vectors.hs new file mode 100644 --- /dev/null +++ b/PDButil/Vectors.hs @@ -0,0 +1,51 @@ +module PDButil.Vectors where + +-- A minimal implementation of vector operations for pdb calculations. +-- These vector operations are NOT safe... they will not verify dimensionality requirements +-- TODO; matrix multiplication, factorizations. + +import Data.List + +dot :: (Num a) => [a] -> [a] -> a +dot a b = foldr1 (+) $ zipWith (*) a b + +-- Only defined on 3 dimensional vectors; no obvious generalization +cross :: (Num a) => [a] -> [a] -> [a] +cross [a1,a2,a3] [b1,b2,b3] = [c1,c2,c3] where + c1 = a2*b3 - a3*b2 + c2 = a3*b1 - a1*b3 + c3 = a1*b2 - a2*b1 + +vAdd :: (Num a) => [a] -> [a] -> [a] +vAdd = zipWith (+) + +vSub :: (Num a) => [a] -> [a] -> [a] +vSub = zipWith (-) + +magnitude :: (Num a) => [a] -> a +magnitude = sum . (map (^2)) + +norm :: (Floating a) => [a] -> a +norm = sqrt . magnitude + +unit :: (Floating a) => [a] -> [a] +unit vec = map (/ (norm vec)) vec + +angle :: (Floating a) => [a] -> [a] -> a +angle a b + | norm a == 0 || norm b == 0 = 0 + | otherwise = acos $ (a `dot` b) / ((norm a) * (norm b)) + +-- Two Dimensions +{- +vRotate :: (Floating a) => a -> [a] -> [a] +vRotate degs [x,y] = mtimes rMatrix [x,y] where + rMatrix = [[cos(degs),-sin(degs)],[sin(degs),cos(degs)]] +-} + +vRotate3d :: (Floating a) => a -> a -> a -> [a] -> [a] +vRotate3d theta phi psi vect = [r1 `dot` vect, r2 `dot` vect, r3 `dot` vect] where + rMatrix = [r1,r2,r3] + r1 = [cos(theta)*cos(psi),-cos(phi)*sin(psi)+sin(phi)*sin(theta)*cos(psi),sin(phi)*sin(psi)+cos(phi)*sin(theta)*cos(psi)] + r2 = [cos(theta)*sin(psi),cos(phi)*cos(psi)+sin(phi)*sin(theta)*sin(psi),-sin(phi)*cos(psi)+cos(phi)*sin(theta)*sin(psi)] + r3 = [-sin(theta),sin(phi)*cos(theta),cos(phi)*cos(theta)] diff --git a/Tests/Test.hs b/Tests/Test.hs --- a/Tests/Test.hs +++ b/Tests/Test.hs @@ -2,8 +2,8 @@ --Try your own PDB as an extra precaution module Main where -import PDBtools.PDButil -import PDBtools.PDBparse +import PDButil.PDBparse +import PDBtools.Base main = do contents <- parse "3C22.pdb"