PDBtools 0.0.2 → 0.0.3
raw patch · 7 files changed
+233/−188 lines, 7 files
Files
- PDBtools.cabal +4/−4
- PDBtools/Base.hs +75/−48
- PDBtools/PDButil/PDBparse.hs +90/−0
- PDBtools/PDButil/Vectors.hs +56/−0
- PDBtools/Residues.hs +8/−2
- PDButil/PDBparse.hs +0/−83
- PDButil/Vectors.hs +0/−51
PDBtools.cabal view
@@ -1,5 +1,5 @@ Name: PDBtools-Version: 0.0.2+Version: 0.0.3 License: GPL-3 License-file: LICENSE Cabal-Version: >= 1.6@@ -16,7 +16,7 @@ Homepage: http://www.github.com/rotskoff Build-type: Simple-Stability: Experiemental+Stability: Experimental extra-source-files: README @@ -33,8 +33,8 @@ exposed-modules: PDBtools.Base PDBtools.Residues- PDButil.PDBparse- PDButil.Vectors+ PDBtools.PDButil.PDBparse+ PDBtools.PDButil.Vectors ghc-options:
PDBtools/Base.hs view
@@ -11,12 +11,11 @@ -- 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 PDBtools.PDButil.PDBparse+import PDBtools.PDButil.Vectors import Data.List import Data.ByteString.Char8 (ByteString) import qualified Data.ByteString.Char8 as B@@ -33,6 +32,10 @@ atomname t atmlist = filter matcht atmlist where matcht a = name a == B.pack t +atomnames :: [String] -> [Atom] -> [Atom]+atomnames ts atmlist = filter matchts atmlist where+ matchts a = name a `elem` (map B.pack ts)+ -- Match the residue type, input the three letter abbreviation restype :: String -> [Atom] -> [Atom] restype t atmlist = filter matcht atmlist where@@ -94,56 +97,80 @@ 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'+dihedrals :: Protein -> [(Double,Double)]+dihedrals p = map (\cAlpha -> dihedral cAlpha (atoms p)) $ backbone p -{---- 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+kth :: Int -> Atom -> String -> [Atom] -> [Double]+kth k cAlpha t atms+ | (getAdjacent k cAlpha t atms) == [] = [0,0,0]+ | otherwise = coords $ head $ getAdjacent k cAlpha t atms --- 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+getAdjacent :: Int -> Atom -> String -> [Atom] -> [Atom]+getAdjacent k cAlpha t atms = filter (\s -> resid s == resid cAlpha+k) $ atomname t atms --} +--TODO TESTING+dihedral :: Atom -> [Atom] -> (Double,Double)+dihedral cAlpha atms+ | prevCO == [0,0,0] || nextN == [0,0,0] = (0,0)+ | name cAlpha == B.pack "CA" = (phi v1 v2,psi v2 v3)+ | otherwise = error "Please input a Carbon Alpha atom." where+ prevCO = kth (-1) cAlpha "C" atms+ currCA = coords $ cAlpha+ currN = kth 0 cAlpha "N" atms+ currCO = kth 0 cAlpha "C" atms+ nextN = kth 1 cAlpha "N" atms+ a = vSub prevCO currCA+ b = vSub currN currCA+ c = vSub currCA currCO+ d = vSub currCA nextN+ v1 = a `cross` b+ v2 = c `cross` b+ v3 = c `cross` d+ phi v1 v2+ | a `dot` v1 < 0 = angle v1 v2+ | otherwise = -(angle v1 v2)+ psi v2 v3+ | d `dot` v3 > 0 = angle v2 v3+ | otherwise = -(angle v2 v3) +rama :: Protein -> IO()+rama p = do+ let toTSV (phi,psi) = (show phi) ++ "\t" ++ (show psi) + mapM_ putStrLn $ map toTSV $ dihedrals p +-- 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) +convert :: ByteString -> String +convert name+ | query == Nothing = "X" + | otherwise = fromJust query where + query = 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'
+ PDBtools/PDButil/PDBparse.hs view
@@ -0,0 +1,90 @@+-- Module : PDBparse+-- Copyright : (c) 2012 Grant Rotskoff+-- License : GPL-3+--+-- Maintainer : gmr1887@gmail.com+-- Stability : experimental++module PDBtools.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)++
+ PDBtools/PDButil/Vectors.hs view
@@ -0,0 +1,56 @@+-- Module : Vectors+-- Copyright : (c) 2012 Grant Rotskoff+-- License : GPL-3+--+-- Maintainer : gmr1887@gmail.com+-- Stability : experimental++module PDBtools.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+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 :: (Num a, Floating a) => [a] -> a+norm = sqrt . magnitude++unit :: (Num a, Floating a) => [a] -> [a]+unit vec = map (/ (norm vec)) vec++angle :: (Eq a, Floating a) => [a] -> [a] -> a+angle a b+ | norm a == 0 || norm b == 0 = 0+ | otherwise = (180/pi) * (acos $ (a `dot` b) / ((norm a) * (norm b)))++rotateAboutz :: (Floating a) => [a] -> a -> [a]+rotateAboutz vec degs = map (dot vec) rMatrix where+ rMatrix = [[cos(theta),sin(theta),0],[-sin(theta),cos(theta),0],[0,0,1]]+ theta = degs*(pi/180)++rotateAboutx :: (Floating a) => [a] -> a -> [a]+rotateAboutx vec degs = map (dot vec) rMatrix where+ rMatrix = [[1,0,0],[0,cos(theta),sin(theta)],[0,-sin(theta),cos(theta)]]+ theta = degs*(pi/180)+++
PDBtools/Residues.hs view
@@ -1,9 +1,15 @@--- Residues centered at the Carbon-Alpha, some sort of directionality constraint+-- Module : Residues+-- Copyright : (c) 2012 Grant Rotskoff+-- License : GPL-3+--+-- Maintainer : gmr1887@gmail.com+-- Stability : experimental + module PDBtools.Residues where import PDBtools.Base-import PDButil.PDBparse+import PDBtools.PDButil.PDBparse --One should really only use these methods on proteins, but for the sake of composing selections, the input form is [Atom]
− PDButil/PDBparse.hs
@@ -1,83 +0,0 @@-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
@@ -1,51 +0,0 @@-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)]