diff --git a/PDBtools.cabal b/PDBtools.cabal --- a/PDBtools.cabal +++ b/PDBtools.cabal @@ -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: diff --git a/PDBtools/Base.hs b/PDBtools/Base.hs --- a/PDBtools/Base.hs +++ b/PDBtools/Base.hs @@ -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' diff --git a/PDBtools/PDButil/PDBparse.hs b/PDBtools/PDButil/PDBparse.hs new file mode 100644 --- /dev/null +++ b/PDBtools/PDButil/PDBparse.hs @@ -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) + + diff --git a/PDBtools/PDButil/Vectors.hs b/PDBtools/PDButil/Vectors.hs new file mode 100644 --- /dev/null +++ b/PDBtools/PDButil/Vectors.hs @@ -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) + + + diff --git a/PDBtools/Residues.hs b/PDBtools/Residues.hs --- a/PDBtools/Residues.hs +++ b/PDBtools/Residues.hs @@ -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] diff --git a/PDButil/PDBparse.hs b/PDButil/PDBparse.hs deleted file mode 100644 --- a/PDButil/PDBparse.hs +++ /dev/null @@ -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) - - diff --git a/PDButil/Vectors.hs b/PDButil/Vectors.hs deleted file mode 100644 --- a/PDButil/Vectors.hs +++ /dev/null @@ -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)]