packages feed

PDBtools 0.0.1 → 0.0.2

raw patch · 8 files changed

+324/−224 lines, 8 files

Files

PDBtools.cabal view
@@ -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:
+ PDBtools/Base.hs view
@@ -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++-}+++++
− PDBtools/PDBparse.hs 
@@ -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)--
− PDBtools/PDButil.hs 
@@ -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---}-----
+ PDBtools/Residues.hs view
@@ -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
+ PDButil/PDBparse.hs view
@@ -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)++
+ PDButil/Vectors.hs view
@@ -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)]
Tests/Test.hs view
@@ -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"