hPDB-examples (empty) → 0.99
raw patch · 13 files changed
+881/−0 lines, 13 filesdep +AC-Vectordep +GLUTdep +Octreesetup-changed
Dependencies added: AC-Vector, GLUT, Octree, OpenGL, QuickCheck, base, bytestring, containers, deepseq, directory, ghc-prim, hPDB, mtl, template-haskell, text, text-format, vector
Files
- LICENSE +31/−0
- Setup.hs +4/−0
- examples/CanonicalAxes.hs +74/−0
- examples/CleanPDB.hs +47/−0
- examples/PDB2Fasta.hs +75/−0
- examples/Rg.hs +51/−0
- examples/ShiftToCenter.hs +30/−0
- examples/SplitModels.hs +34/−0
- examples/StericClashCheck.hs +49/−0
- examples/Viewer.hs +275/−0
- hPDB-examples.cabal +87/−0
- tests/PrintEvents.hs +113/−0
- tests/PrintStructureObject.hs +11/−0
+ LICENSE view
@@ -0,0 +1,31 @@+PDB data in here (*.pdb files) is subject to Protein Databank License.++Haskell code in this package is subject to:++Copyright (c) Michal J. Gajda 2010-2013++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions+are met:+1. Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.+2. Redistributions in binary form must reproduce the above copyright+ notice, this list of conditions and the following disclaimer in the+ documentation and/or other materials provided with the distribution.+3. Neither the name of the author nor the names of his contributors+ may be used to endorse or promote products derived from this software+ without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE+ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS+OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)+HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT+LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY+OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF+SUCH DAMAGE.
+ Setup.hs view
@@ -0,0 +1,4 @@+#! /usr/bin/env runhaskell++import Distribution.Simple+main = defaultMain
+ examples/CanonicalAxes.hs view
@@ -0,0 +1,74 @@+{-# LANGUAGE OverloadedStrings, BangPatterns, FlexibleInstances,UndecidableInstances #-}++module Main where++import Control.Monad(when)+import System.IO+import System.Environment+import System.Exit(exitFailure, exitSuccess)+import qualified Bio.PDB.IO as PDBIO+import qualified Data.ByteString.Char8 as BS+import Bio.PDB.Iterable as It+import Bio.PDB.Structure as PDBS+import Bio.PDB.Structure.Vector+import Data.List+import Text.Printf++ifoldrPairs fred fpair e s = pairs+ where+ pairs' a cont = It.ifoldr (\at r -> (at `fpair` a) `fred` r) cont (s :: Structure)+ pairs = It.ifoldr (\at r -> pairs' at r) e (s :: Structure)++ifoldPairs fred fpair e s = pairs+ where+ pairs' a cont = It.ifoldl' (\r at -> (at `fpair` a) `fred` r) cont (s :: Structure)+ pairs = It.ifoldl' (\r at -> pairs' at r) e (s :: Structure)++-- | findAxes finds all three principal axes so that dimensions are ordered.+-- TODO: find a way to loop v1, v2, v3 assignment+findAxes structure = let v1 = findLongestOrthogonalVector [ ] structure+ axis1 = vnormalise v1+ dim1 = vnorm v1+ v2 = findLongestOrthogonalVector [axis1 ] structure+ axis2 = vnormalise v2+ dim2 = vnorm v2+ v3 = findLongestOrthogonalVector [axis1, axis2] structure+ axis3 = vnormalise v3+ dim3 = vnorm v3+ in dim1 `seq` dim2 `seq` dim3 `seq` ([dim1, dim2, dim3], [axis1, axis2, axis3])+ where+ findLongestOrthogonalVector axes = ifoldPairs pickMaxDist (atDistPerpend axes) nullVector + nullVector = fromInteger 0+ atDistPerpend axes !a1 !a2 = vperpends ((coord a1) - (coord a2)) axes+ pickMaxDist !v1 !v2 = if vnorm v1 > vnorm v2 then v1 else v2+ ++{-center :: Structure -> (Double, Double, Double)+center s = (realToFrac sx, realToFrac sy, realToFrac sz)+ where+ + result = It.ifoldr + maxFst (a, b) (c, d) = if a >= c then (a, b) else (c, d)+ coords = [It.ifoldr (\at r -> coord (at :: PDBS.Atom) : r) [] (s :: Structure)]+-}+++main = do args <- getArgs+ when (length args /= 2) $ do hPutStrLn stderr "USAGE: CanonicalAxes <input.pdb> <output.pdb>"+ exitFailure+ let [inpfname, outfname] = args+ Just structure <- PDBIO.parse $ BS.pack inpfname+ let ([d1, d2, d3], axes@[yaxis, xaxis, zaxis]) = findAxes structure+ printf "Dimensions: %.2f %.2f %.2f\n" d1 d2 d3+ putStr "Axis 1 (Y): "+ print yaxis+ putStr "Axis 2 (X): "+ print xaxis+ putStr "Axis 3 (Z): "+ print zaxis+ {- This requires implementation of geometric transforms+ let xform = axesToTransform axes+ let s2 = applyTransform xform s1+ PDBIO.write s2 $ BS.pack outfname+ -}+ exitSuccess
+ examples/CleanPDB.hs view
@@ -0,0 +1,47 @@+{-# LANGUAGE OverloadedStrings, FlexibleInstances, UndecidableInstances, ScopedTypeVariables, FlexibleContexts #-}++module Main where++import System.Environment+import Bio.PDB as PDB+import qualified Data.ByteString.Char8 as BS+import Data.List+import Text.Printf+import Data.Vector.V3+import Control.Monad.State(State, modify, get, runState)++type CounterM a = State Int a++counter :: CounterM Int+counter = do r <- get+ modify (+1)+ return r++runCounterM :: CounterM a -> a+runCounterM = fst . (flip runState) 1++renumberResidues :: (Iterable s Chain) => s -> s+renumberResidues = imap (\ch -> runCounterM .+ imapM forEachResidue $ (ch :: Chain))+ where+ forEachResidue r = do v <- counter+ return $ r { resSeq = v }++renumberAtoms :: (Iterable s Model) => s -> s+renumberAtoms = imap (\m -> runCounterM $+ imapM forEachChain (m :: Model))+ where+ forEachChain :: Chain -> CounterM Chain+ forEachChain ch = do newCh <- imapM forEachAtom ch+ counter+ return newCh+ forEachAtom :: Atom -> CounterM Atom+ forEachAtom at = do v <- counter+ return $ at { atSerial = v }++main = do [inpfname, outfname] <- getArgs+ Just structure <- PDB.parse $ BS.pack inpfname+ putStrLn $ show (PDB.numAtoms structure) ++ " atoms."+ let s1 = renumberResidues $ renumberAtoms $ structure+ PDB.write s1 $ BS.pack outfname+
+ examples/PDB2Fasta.hs view
@@ -0,0 +1,75 @@+{-# LANGUAGE OverloadedStrings, FlexibleContexts, NoMonomorphismRestriction #-}+module Main(main) where++import System.Environment(getProgName, getArgs)+import System.Exit(exitWith, ExitCode(ExitSuccess), ExitCode(..))+import System.Console.GetOpt+import Bio.PDB as PDB+import Bio.PDB.Fasta(fastaRecord)+import qualified Data.ByteString.Char8 as BS+import qualified Data.Map as Map+import Bio.PDB.Fasta+import Control.Monad(when)++data Options = Options { gapped :: Bool,+ allModels :: Bool+ }++showHelp :: IO ()+showHelp =+ do prg <- getProgName+ Prelude.putStrLn $ usageInfo prg options++showVersion = Prelude.putStrLn "Version 0.1" -- TODO: extract version from CABAL declaration++exitAfter function exitCode _ =+ do function+ exitWith exitCode+ return undefined++defaultOptions = Options { gapped = False+ , allModels = False+ }++options :: [OptDescr (Options -> IO Options)]+options = + [Option "V" ["version"] (NoArg (exitAfter showVersion ExitSuccess ))+ "Print program version.",+ Option "h" ["help"] (NoArg (exitAfter showHelp ExitSuccess ))+ "Prints help.",+ Option "g" ["with-gaps"] (NoArg (\opts -> return $ opts { gapped = True }))+ "Shows sequence with gaps.",+ Option "a" ["all-models"] (NoArg (\opts -> return $ opts { allModels = True }))+ "Shows sequences for all MODELs."+ ]++main = do+ args <- getArgs+ let (actions, filenames, opts) = getOpt RequireOrder options args+ opts <- foldl (>>=) (return defaultOptions) actions+ when (length filenames == 0) $ exitAfter showHelp (ExitFailure 3) undefined+ mapM_ (processFile opts) filenames+ exitWith ExitSuccess -- TODO: only if no errors!+-- | This program loads a PDB file and outputs FASTA sequences of all chains within++printFastaRecords :: Options -> [Char] -> Structure -> IO ()+printFastaRecords opts fname s = if allModels opts+ then+ process s+ else+ case firstModel s of+ Nothing -> Prelude.putStrLn $ "No models found within file '" ++ fname ++ "'."+ Just m -> process m+ where+ process = PDB.ifoldM (\() ch -> Prelude.putStrLn $ mkRecord fname ch) ()+ mkRecord = if gapped opts+ then fastaGappedRecord+ else fastaRecord ++--printFastaRecords fname s = PDB.ifoldM (\() ch -> Prelude.putStrLn $ fastaRecord fname ch) () s++processFile opts fname = do maybePDB <- PDB.parse $ BS.pack fname+ case maybePDB of+ Nothing -> return ()+ Just structure -> printFastaRecords opts fname structure+
+ examples/Rg.hs view
@@ -0,0 +1,51 @@+{-# LANGUAGE ScopedTypeVariables, OverloadedStrings, BangPatterns, FlexibleContexts, OverloadedStrings #-}+module Main(main) where++import Bio.PDB as PDB+import System.Environment(getArgs)+import System.IO.Unsafe(unsafePerformIO)+import System.IO(stderr, stdout, hFlush)+import Data.ByteString.Char8 as BS+import Control.DeepSeq+import Text.Printf+import System.Exit+--import Control.Parallel++import Bio.PDB.Structure.Elements(atomicMass, assignElement)++main = do args <- getArgs+ mapM (readAndComputeRg . BS.pack) args+ exitWith ExitSuccess++readAndComputeRg filename =+ do Just structure <- PDB.parse filename+ BS.putStr "Parsed "+ hFlush stdout+ BS.putStr . BS.pack . show $ numAtoms structure+ hFlush stdout+ rnf structure `seq` BS.putStrLn " atoms!"+ rg <- return $! radiusOfGyration structure+ printf "%s: %.2f\n" (BS.unpack filename) rg + return rg++-- TODO: consider mass of an atom+-- TODO: use Data.Tensor or some other standard notation for 3D vectors+-- TODO: check if there is nice library with quaternions and rotation matrices for these+radiusOfGyration :: Iterable Structure Atom => Structure -> Double+radiusOfGyration structure = avgDistDev+ where+ -- (c -> b -> c) -> c -> a -> c+ avgDistDev = sqrt ((ifoldl' addDistDev 0.0 structure)/totalMass)+ addDistDev !total at = total + (vnorm $ coord at - center)**2 * atMass at+ atMass :: Atom -> Double+ atMass at = atomicMass $ case assignElement at of+ "" -> "C"+ otherwise -> otherwise+ center = v |* (1.0/totalMass)+ where v = ifoldl' addCoord nullVector structure+ nullVector = fromInteger 0+ addCoord v (at :: Atom) = v + coord at |* atMass at+ counter x (at :: Atom) = x + atMass at + totalMass = ifoldl' counter 0.0 structure++
+ examples/ShiftToCenter.hs view
@@ -0,0 +1,30 @@+{-# LANGUAGE OverloadedStrings, FlexibleInstances, UndecidableInstances #-}++module Main where++import System.Environment+import Bio.PDB as PDB+import qualified Data.ByteString.Char8 as BS+import Data.List+import Text.Printf+import Data.Vector.V3++center :: PDB.Structure -> Vector3+center s = avgv+ where+ sumv = ifoldl' addCoord 0 (s :: PDB.Structure)+ n = realToFrac . fromIntegral . PDB.numAtoms $ s+ addCoord v (PDB.Atom { coord = c }) = v+c+ avgv = (1/n) *| sumv++shift v = PDB.imap subCoord+ where+ subCoord (at@Atom { coord = c }) = at { coord = c - v }++main = do [inpfname, outfname] <- getArgs+ Just structure <- PDB.parse $ BS.pack inpfname+ let c@(Vector3 x y z) = center structure+ printf "Center %.2f %.2f %.2f\n" x y z+ putStrLn $ show (PDB.numAtoms structure) ++ " atoms."+ let s1 = shift c structure+ PDB.write s1 $ BS.pack outfname
+ examples/SplitModels.hs view
@@ -0,0 +1,34 @@+{-# LANGUAGE OverloadedStrings, FlexibleInstances, UndecidableInstances #-}++module Main where++import System.Environment+import System.Exit+import System.IO(hPutStrLn, stderr)+--import qualified Bio.PDB.IO as PDBIO+import qualified Data.ByteString.Char8 as BS+--import Bio.PDB.Iterable as It+import Bio.PDB as PDB+import Data.List+import Text.Printf+--import Data.Vector.V3+import qualified Data.Vector as V+import Control.Monad(forM, when)+--import Control.DeepSeq++splitModels aStructure = map mkStructure . V.toList . models $ aStructure+ where mkStructure aModel = aStructure{ models = V.singleton aModel }++usage = do hPutStrLn stderr $ concat [ "Usage: SplitModels <input.pdb> <output_prefix>"+ , "Splitting of PDB files with multiple MODEL entries." ]+ exitFailure++main = do lenArgs <- length `fmap` getArgs+ when (lenArgs /= 2) usage+ [inpfname, outfname] <- getArgs+ Just structure <- PDB.parse $ BS.pack inpfname+ let splitted = zip [1..] $ splitModels structure+ forM splitted $ \(num, aStructure) ->+ do let fname = outfname ++ show num ++ ".pdb"+ putStrLn $ "Writing " ++ show fname ++ " with " ++ show (numAtoms aStructure) ++ " atoms."+ PDB.write aStructure $ BS.pack fname
+ examples/StericClashCheck.hs view
@@ -0,0 +1,49 @@+{-# LANGUAGE OverloadedStrings, DisambiguateRecordFields, ScopedTypeVariables #-}++module Main where++import qualified System.Environment as Env+import qualified Data.ByteString.Char8 as BS+import qualified Data.Octree as Oct+import Bio.PDB as PDB+import qualified Bio.PDB.Structure.Elements as Elt++-- TODO: way to hand over Atom object along with its Residue, Chain and Model+-- TODO: make it so that Show gives a "full_id" string by default, e.g. Model, chain, residue id, atom name.+clashCheck s1 s2 = filter (/= []) . Prelude.map clashes $ ifoldr (:) [] s2+ where+ clashes (at :: PDB.Atom) = Oct.withinRange ot (radius + maxRadius) (v3v $ PDB.coord at)+ where+ radius :: Double = realToFrac $ Elt.vanDerWaalsRadius . PDB.element $ at+ ot :: Oct.Octree (Int, Double)+ ot = makeOctree s1++extract :: PDB.Atom -> (Oct.Vector3, (Int, Double))+extract (PDB.Atom { coord = cvec,+ atSerial = ser ,+ element = elt }) = (v3v cvec,+ (ser, realToFrac vdw))+ where+ vdw = Elt.vanDerWaalsRadius elt++makeOctree structure = Oct.fromList . Prelude.map extract . ifoldr (:) [] $ structure++-- Bio.PDB.Structure.Elements should export max bound for vdw etc.+-- or a list of known element codes.+maxRadius = 1.6++v3v = id++size ot = length . Oct.toList $ ot++main = do [input1, input2] <- Env.getArgs+ Just structure1 <- PDB.parse $ BS.pack input1+ Just structure2 <- PDB.parse $ BS.pack input2+ print $ clashCheck structure1 structure2+ print "Depths:"+ print . Oct.depth . makeOctree $ structure1+ print . Oct.depth . makeOctree $ structure2+ print "Sizes:"+ print . size . makeOctree $ structure1+ print . size . makeOctree $ structure2+
+ examples/Viewer.hs view
@@ -0,0 +1,275 @@+{-# LANGUAGE ScopedTypeVariables, BangPatterns #-}+import System.IO.Unsafe(unsafePerformIO)+import Data.IORef as IORef+import Control.Monad(ap)++-- OpenGL stuff +import qualified Graphics.Rendering.OpenGL as GL+import qualified Graphics.UI.GLUT as GL+import Graphics.Rendering.OpenGL(($=))++-- PDB reading stuff+import Bio.PDB.Structure as PDBS+import Bio.PDB.Structure.Elements as PDBE+import Bio.PDB.IO as PDBIO+import Bio.PDB.Iterable as PDBI+import Data.ByteString.Char8 as BS+import Bio.PDB.Structure.Vector++data UIState = UIState { manipulationCenter :: Maybe GL.Position+ , pressedButton :: Maybe GL.MouseButton+ , manipulatedMatrix :: Maybe (GL.GLmatrix GL.GLdouble)+ , windowSize :: GL.Position+ , currentViewMatrix :: GL.GLmatrix GL.GLdouble --[GL.GLdouble]+ }++initialUIState = do GL.loadIdentity+ initialViewMatrix :: GL.GLmatrix GL.GLdouble <- GL.get $ GL.matrix $ Nothing+ --m <- GL.getMatrixComponents GL.ColumnMajor initialViewMatrix+ return $ cleanUIState initialViewMatrix++cleanUIState initialViewMatrix = UIState { manipulationCenter = Nothing+ , pressedButton = Nothing+ , manipulatedMatrix = Nothing+ , windowSize = GL.Position 1 1+ , currentViewMatrix = initialViewMatrix+ }++setWindowSize :: GL.GLsizei -> GL.GLsizei -> UIState -> UIState+setWindowSize w h uiState = uiState { windowSize = GL.Position (fromIntegral w) (fromIntegral h) }++main :: IO ()+main = do+ (progname, [filename]) <- GL.getArgsAndInitialize+ -- First read structure+ Just structure <- PDBIO.parse $ BS.pack filename+ -- Initialize OpenGL+ --displayMode $= [With DisplayDepth, With DisplayRGB] -- is it glutInitDisplayModel (GLUT_SINGLE | GLUT_RGB | GLUT_DEPTH)+ setup progname -- setup OpenGL scene model+ -- Initialize camera+ initUIState <- initialUIState+ uiState <- newIORef initUIState+ translateViewMatrix uiState $ vec3DToVector3 $ (-(center structure))+ translateViewMatrix uiState $ GL.Vector3 0 0 (-frustumMiddle)+ -- Now set callbacks+ GL.displayCallback $= display uiState structure+ GL.reshapeCallback $= Just (reshape uiState)+ GL.keyboardMouseCallback $= Just (handleKeys uiState)+ GL.motionCallback $= Just (handleMouseMotion uiState)+ GL.closeCallback $= Just (appClosing uiState)+ display uiState structure+ GL.mainLoop++appClosing uiState = do uiSt <- readIORef uiState+ let mat = currentViewMatrix uiSt+ cm <- GL.getMatrixComponents GL.ColumnMajor mat+ Prelude.putStr "Final view matrix:"+ Prelude.print cm+++modifyViewMatrix :: IORef UIState -> IO () -> IO ()+modifyViewMatrix uiState modification =+ do uiSt <- GL.get uiState+ GL.matrixMode $= GL.Modelview 0+ GL.matrix Nothing $= currentViewMatrix uiSt+ --m :: GL.GLmatrix GL.GLdouble <- GL.newMatrix GL.ColumnMajor $ currentViewMatrix uiSt+ --GL.matrix Nothing $= m+ modification+ m :: GL.GLmatrix GL.GLdouble <- GL.get $ GL.matrix Nothing+ --cm <- GL.getMatrixComponents GL.ColumnMajor m+ --print cm+ uiState $= uiSt { currentViewMatrix = m }+ GL.postRedisplay Nothing++manipulateViewMatrix :: IORef UIState -> IO () -> IO (GL.GLmatrix GL.GLdouble)+manipulateViewMatrix uiState modification =+ do uiSt <- GL.get uiState+ GL.matrixMode $= GL.Modelview 0+ GL.matrix Nothing $= unJust (manipulatedMatrix uiSt)+ --GL.matrix Nothing $= unJust (manipulatedMatrix uiSt)+ --GL.loadIdentity+ --GL.multMatrix $ unJust $ manipulatedMatrix uiSt+ modification+ --GL.matrix Nothing $= unJust (manipulatedMatrix uiSt)+ m :: GL.GLmatrix GL.GLdouble <- GL.get $ GL.matrix Nothing+ cm <- GL.getMatrixComponents GL.ColumnMajor m+ --print cm+ uiState $= uiSt { currentViewMatrix = m }+ GL.postRedisplay Nothing+ return m++xaxis = GL.Vector3 1 0 0 :: GL.Vector3 GL.GLdouble++yaxis = GL.Vector3 0 1 0 :: GL.Vector3 GL.GLdouble++zaxis = GL.Vector3 0 0 1 :: GL.Vector3 GL.GLdouble++negateVector3 :: (Num a) => GL.Vector3 a -> GL.Vector3 a+negateVector3 = fmap (\a -> (-a))++handleKeys :: IORef UIState -> GL.Key -> GL.KeyState -> GL.Modifiers -> GL.Position -> IO ()+handleKeys uiState (GL.Char 'q') GL.Down _modifiers _position = rotateViewMatrix uiState 10.0 yaxis+handleKeys uiState (GL.Char 'e') GL.Down _modifiers _position = rotateViewMatrix uiState (-10.0) yaxis+handleKeys uiState (GL.Char 'r') GL.Down _modifiers _position = rotateViewMatrix uiState 10.0 xaxis+handleKeys uiState (GL.Char 'f') GL.Down _modifiers _position = rotateViewMatrix uiState (-10.0) xaxis+handleKeys uiState (GL.Char 'a') GL.Down _modifiers _position = translateViewMatrix uiState xaxis+handleKeys uiState (GL.Char 'd') GL.Down _modifiers _position = translateViewMatrix uiState $ negateVector3 xaxis+handleKeys uiState (GL.Char 'w') GL.Down _modifiers _position = translateViewMatrix uiState zaxis+handleKeys uiState (GL.Char 's') GL.Down _modifiers _position = translateViewMatrix uiState $ negateVector3 zaxis+handleKeys uiState (GL.MouseButton GL.RightButton) GL.Down _modifiers position = startMouseTracking uiState GL.RightButton position +handleKeys uiState (GL.MouseButton GL.MiddleButton) GL.Down _modifiers position = startMouseTracking uiState GL.MiddleButton position +handleKeys uiState (GL.MouseButton GL.LeftButton) GL.Down _modifiers position = startMouseTracking uiState GL.LeftButton position +handleKeys uiState (GL.MouseButton GL.MiddleButton) GL.Up _modifiers position = do matrix <- mouseMoveViewMatrix xyAxes uiState position+ modifyIORef uiState $ clearManipulatedMatrix matrix+handleKeys uiState (GL.MouseButton GL.RightButton ) GL.Up _modifiers position = do matrix <- mouseMoveViewMatrix xzAxes uiState position+ modifyIORef uiState $ clearManipulatedMatrix matrix+handleKeys uiState (GL.MouseButton GL.LeftButton) GL.Up _modifiers position = do matrix <- mouseRotateViewMatrix uiState position+ modifyIORef uiState $ clearManipulatedMatrix matrix+handleKeys _ _ _ _ _ = return ()++startMouseTracking uiState button position = modifyIORef uiState newState+ where+ newState uiState = uiState { manipulationCenter = Just position+ , pressedButton = Just button+ , manipulatedMatrix = Just $ currentViewMatrix uiState+ }+++unJust (Just a) = a+unJust Nothing = error "unJust of Nothing"++rotateViewMatrix :: IORef UIState -> GL.GLdouble -> GL.Vector3 GL.GLdouble -> IO ()+rotateViewMatrix uiState angle axis =+ modifyViewMatrix uiState $ GL.rotate angle axis++clearManipulatedMatrix matrix uiState = uiState { manipulationCenter = Nothing+ , pressedButton = Nothing+ , manipulatedMatrix = Nothing+ , currentViewMatrix = matrix+ }++translateViewMatrix :: IORef UIState -> GL.Vector3 GL.GLdouble -> IO ()+translateViewMatrix uiState translation =+ modifyViewMatrix uiState $ GL.translate translation++mouseMoveViewMatrix axesExpansion uiState pos1 = do uiSt <- readIORef uiState+ let pos2 = unJust $ manipulationCenter uiSt+ let windowDim = windowSize uiSt+ let (x, y) = getPosChange pos1 pos2 windowDim+ manipulateViewMatrix uiState $ GL.translate $ axesExpansion x y++handleMouseMotion uiState position = do btn <- pressedButton `fmap` readIORef uiState+ case btn of+ Just GL.MiddleButton -> do mouseMoveViewMatrix xyAxes uiState position+ return ()+ Just GL.RightButton -> do mouseMoveViewMatrix xzAxes uiState position+ return ()+ Just GL.LeftButton -> do mouseRotateViewMatrix uiState position+ return ()++xyAxes x y = GL.Vector3 (x * frustumMiddle) (y * frustumMiddle) 0+xzAxes x y = GL.Vector3 (x * frustumMiddle) 0 ((-y) * frustumMiddle)++mouseRotateViewMatrix :: IORef UIState -> GL.Position -> IO (GL.GLmatrix GL.GLdouble)+mouseRotateViewMatrix uiState pos1 = do uiSt <- readIORef uiState+ let pos2 = unJust $ manipulationCenter uiSt+ let windowDim = windowSize uiSt+ let (angle, axis) = computeRotation pos1 pos2 windowDim+ --Prelude.putStrLn $ "AA" ++ show (angle, axis)+ manipulateViewMatrix uiState $ --do --GL.translate $ GL.Vector3 0 0 (-frustumMiddle)+ GL.rotate angle axis+ --GL.translate $ GL.Vector3 0 0 ( frustumMiddle)++computeRotation :: GL.Position -> GL.Position -> GL.Position -> (GL.GLdouble, GL.Vector3 GL.GLdouble)+computeRotation (GL.Position curX curY )+ (GL.Position initX initY)+ (GL.Position w h ) = (angle, GL.Vector3 dy dx 0.0)+ where+ [xf, yf, wf, hf] = Prelude.map fromIntegral [curX - initX, curY - initY, w, h]+ dx = xf / (min hf wf)+ dy = yf / (min hf wf)+ angle = 180.0 * sqrt (dx*dx + dy*dy)+ +getPosChange :: GL.Position -> GL.Position -> GL.Position -> (GL.GLdouble, GL.GLdouble)+getPosChange (GL.Position x y) (GL.Position x' y') (GL.Position w h) = (fromIntegral (x - x')/fromIntegral w,+ fromIntegral (y' - y )/fromIntegral h)++reflector = GL.Light 0++setup :: Prelude.String -> IO ()+setup progname = do+ GL.initialDisplayMode $= [GL.RGBMode, GL.WithDepthBuffer, GL.DoubleBuffered]+ GL.createWindow progname+ GL.clearColor $= GL.Color4 0 0 0 0+ GL.shadeModel $= GL.Smooth+ GL.materialSpecular GL.FrontAndBack $= GL.Color4 1 1 1 1+ GL.materialShininess GL.FrontAndBack $= 50+ GL.position reflector $= GL.Vertex4 1 1 1 0+ GL.lighting $= GL.Enabled+ GL.light reflector $= GL.Enabled+ GL.depthFunc $= Just GL.Less+ +-- Assessing dimensions for initial focus+center :: PDBS.Structure -> Vec3D+center structure = average + where+ (!average, _count) = PDBI.ifoldl' step (fromIntegral 0, 0) structure+ step :: (Vec3D, Double) -> PDBS.Atom -> (Vec3D, Double)+ step (!r, !i) at = let i' = i + 1+ in (coord at |* (1/i') + r |* (i/i'), i')++dims structure = maxv - minv+ where+ (!minv, !maxv) = PDBI.ifoldl' (\(!minv, !maxv) at -> let c = coord (at :: PDBS.Atom)+ in (vzip min minv c,+ vzip max maxv c)) (cs, cs) structure+ !cs = center structure++vec3DToVector3 :: Vec3D -> GL.Vector3 GL.GLdouble+vec3DToVector3 v = GL.Vector3 x' y' z'+ where+ (x, y, z) = unpackVec3D v+ [x', y', z'] :: [GL.GLdouble] = Prelude.map realToFrac [x, y, z]++renderStructure :: IORef UIState -> PDBS.Structure -> IO ()+renderStructure uiState structure = PDBI.ifoldM (\_ -> renderAtom) () structure+ where renderAtom (at :: Atom) = GL.preservingMatrix $ do GL.matrixMode $= GL.Modelview 0+ uiSt <- readIORef uiState+ --m :: GL.GLmatrix GL.GLdouble <- GL.newMatrix GL.ColumnMajor $ currentViewMatrix uiSt+ GL.matrix Nothing $= currentViewMatrix uiSt+ --cm <- GL.getMatrixComponents GL.ColumnMajor $ currentViewMatrix uiSt+ --print cm+ GL.translate $ vec3DToVector3 $ PDBS.coord at+ GL.renderObject GL.Solid $ GL.Sphere' (realToFrac radius) 40 32 -- much smoother looks!+ where+ radius = realToFrac $ PDBE.vanDerWaalsRadius $ PDBS.element at++display :: IORef UIState -> PDBS.Structure -> IO ()+display uiState structure = do GL.clear [GL.ColorBuffer,+ GL.DepthBuffer]+ renderStructure uiState structure+ GL.swapBuffers++myOrtho w h = let hw = fromIntegral h/fromIntegral w+ wh = fromIntegral w/fromIntegral h+ in if w <= h+ then GL.ortho (-1.5) (1.5) (-1.5*hw) (1.5*hw) (-10) (10)+ else GL.ortho (-1.5*wh) (1.5*wh) (-1.5) (1.5) (-10) (10)++frustumFront = 5.0+frustumBack = 1000.0+frustumMiddle = frustumFront + abs (frustumBack - frustumFront)/2.0+--frustumMiddle = abs (frustumFront - frustumBack) / 2.0+++reshape uiState s@(GL.Size w h) = do GL.viewport $= (GL.Position 0 0, s)+ modifyIORef uiState $ setWindowSize w h+ GL.matrixMode $= GL.Projection+ GL.loadIdentity+ GL.frustum (-x) x (-y) y frustumFront frustumBack+ GL.matrixMode $= GL.Modelview 0+ where+ x = fromIntegral w/fromIntegral h+ y = 1.0++
+ hPDB-examples.cabal view
@@ -0,0 +1,87 @@+name: hPDB-examples+version: 0.99+stability: beta+homepage: https://github.com/mgajda/hpdb+package-url: http://hackage.haskell.org/package/hPDB+synopsis: Parser, print and manipulate structures in PDB file format.+description: Protein Data Bank file format is a most popular format for holding biomolecule data. This is a very fast parser (below 7s for the largest entry in PDB - 1HTQ which is over 70MB - as compared with 11s of RASMOL 2.7.5, or 2m15s of BioPython with Python 2.6 interpreter.) It is aimed to not only deliver event-based interface, but also a high-level data structure for manipulating data in spirit of BioPython's PDB parser. +category: Bioinformatics +license: BSD3+license-file: LICENSE++author: Michal J. Gajda+copyright: Copyright by Michal J. Gajda '2009-'2012+maintainer: mjgajda@googlemail.com+bug-reports: mailto:mjgajda@googlemail.com++build-type: Simple+cabal-version: >=1.8+tested-with: GHC==7.4.1, GHC==7.0.3, GHC==7.4.2+--Need to re-test: GHC==6.12.1, GHC==7.0.4, GHC==7.1.20101026 +--data-files: 1CRN.pdb, 1HTQ.pdb, 3JYV.pdb++source-repository head+ type: git+ location: git://github.com:mgajda/hpdb-examples.git++Executable PDB2Fasta+ main-is: examples/PDB2Fasta.hs+ ghc-options: -fspec-constr-count=4 -O3 + build-depends: base>=4.0, base <4.7, bytestring, ghc-prim, directory, mtl, template-haskell, vector, AC-Vector, containers, deepseq, QuickCheck >= 2.5.0.0, text>=0.11.1.13, hPDB >= 0.99+ other-extensions: ScopedTypeVariables OverloadedStrings BangPatterns NoMonomorphismRestriction EmptyDataDecls MagicHash++Executable ShiftToCenter+ main-is: examples/ShiftToCenter.hs+ ghc-options: -fspec-constr-count=4 -O3 + build-depends: base>=4.0, base <4.7, ghc-prim, directory, mtl, template-haskell, vector, AC-Vector, containers, deepseq,QuickCheck >= 2.5.0.0, text>=0.11.1.13, hPDB >= 0.99, bytestring+ other-extensions: ScopedTypeVariables OverloadedStrings BangPatterns NoMonomorphismRestriction EmptyDataDecls MagicHash++Executable CleanPDB+ main-is: examples/CleanPDB.hs+ ghc-options: -fspec-constr-count=4 -O3 + build-depends: base>=4.0, base <4.7, bytestring, ghc-prim, directory, mtl, template-haskell, vector, AC-Vector, containers, deepseq,QuickCheck >= 2.5.0.0, text>=0.11.1.13, hPDB >= 0.99+ other-extensions: ScopedTypeVariables OverloadedStrings BangPatterns NoMonomorphismRestriction EmptyDataDecls MagicHash++Executable SplitModels+ main-is: examples/SplitModels.hs+ ghc-options: -fspec-constr-count=4 -O3 + build-depends: base>=4.0, base <4.7, ghc-prim, directory, mtl, template-haskell, vector, AC-Vector, containers, deepseq,QuickCheck >= 2.5.0.0, text>=0.11.1.13, hPDB >= 0.99, bytestring+ other-extensions: ScopedTypeVariables OverloadedStrings BangPatterns NoMonomorphismRestriction EmptyDataDecls MagicHash++Executable CanonicalAxes+ main-is: examples/CanonicalAxes.hs+ ghc-options: -fspec-constr-count=4 -O3 + Build-depends: base>=4.0, base <4.7, bytestring, ghc-prim, directory, mtl, template-haskell, vector, AC-Vector, containers, deepseq,QuickCheck >= 2.5.0.0, text>=0.11.1.13, hPDB >= 0.99+ other-extensions: ScopedTypeVariables OverloadedStrings BangPatterns NoMonomorphismRestriction EmptyDataDecls MagicHash++Executable PrintEvents+ Main-is: tests/PrintEvents.hs+ ghc-options: -fspec-constr-count=4 -O3 + Build-depends: base>=4.0, base <4.7, bytestring, ghc-prim, directory, mtl, template-haskell, vector, AC-Vector, containers, deepseq, QuickCheck >= 2.5.0.0, text>=0.11.1.13, hPDB >= 0.99+ other-extensions: ScopedTypeVariables OverloadedStrings BangPatterns NoMonomorphismRestriction EmptyDataDecls++Executable PrintStructureObject+ ghc-options: -fspec-constr-count=4 -O3 + Build-depends: base>=4.0, base <4.7, bytestring, ghc-prim, directory, mtl, template-haskell, vector, AC-Vector, containers, deepseq,QuickCheck >= 2.5.0.0, text>=0.11.1.13, hPDB >= 0.99+ Main-is: tests/PrintStructureObject.hs+ other-extensions: ScopedTypeVariables OverloadedStrings BangPatterns NoMonomorphismRestriction EmptyDataDecls++Executable Rg+ main-is: examples/Rg.hs+ ghc-options: -fspec-constr-count=4 -O3 + Build-depends: base>=4.0, base <4.7, ghc-prim, bytestring, directory, mtl, template-haskell, vector, AC-Vector, containers, deepseq,QuickCheck >= 2.5.0.0, text>=0.11.1.13, hPDB >= 0.99+ other-extensions: ScopedTypeVariables OverloadedStrings BangPatterns NoMonomorphismRestriction EmptyDataDecls++Executable StericClashCheck+ main-is: examples/StericClashCheck.hs+ ghc-options: -fspec-constr-count=4 -O3 + Build-Depends: base>=4.0, base<4.7, bytestring, ghc-prim, directory, mtl, template-haskell, vector, AC-Vector, containers, deepseq, Octree >= 0.5, QuickCheck, text, text-format, hPDB >= 0.4+ other-extensions: ScopedTypeVariables OverloadedStrings BangPatterns NoMonomorphismRestriction EmptyDataDecls++-- if flag(haveOpenGL) -- why can't we make it conditionally available?+Executable Viewer+ main-is: examples/Viewer.hs+ ghc-options: -fspec-constr-count=4 -O3 + Build-depends: base>=4.0, base <4.7, ghc-prim, directory, mtl, template-haskell, vector, AC-Vector, containers, deepseq,QuickCheck >= 2.5.0.0, text>=0.11.1.13, OpenGL, GLUT, hPDB >= 0.4, bytestring+ other-extensions: ScopedTypeVariables OverloadedStrings BangPatterns NoMonomorphismRestriction EmptyDataDecls MagicHash+
+ tests/PrintEvents.hs view
@@ -0,0 +1,113 @@+{-# LANGUAGE OverloadedStrings #-}++module Main where++import System.Environment+import System.Exit+import System.Console.GetOpt+import System.Mem(performGC)+import qualified Control.Exception(catch)+import Control.Exception.Base(SomeException)+import System.IO+--import Bio.PDB.EventParser.PDBParsingAbstractions+import Control.Monad(when)++import qualified Bio.PDB.EventParser.PDBEventPrinter as PDBEventPrinter+import qualified Bio.PDB.EventParser.PDBEvents as PDBEvents+import Bio.PDB.EventParser.PDBEventParser(parsePDBRecords)++import qualified Data.ByteString.Char8 as BS+import Bio.PDB.IO.OpenAnyFile as OpenAnyFile++data Options = Options { optVerbosity :: Int,+ eventFilter :: String -> PDBEvents.PDBEvent -> Maybe String,+ printing :: Bool+ }++showHelp :: IO ()+showHelp =+ do prg <- getProgName+ putStrLn $ usageInfo prg options++showVersion = putStrLn "Version 0.1"++exitAfter function exitCode opts =+ do function+ exitWith exitCode+ return opts++defaultOptions = Options { optVerbosity = 0,+ eventFilter = \f a -> Just $ show a,+ printing = False+ }++changeVerbosity :: (Int -> Int) -> Options -> IO Options+changeVerbosity aChange opt = return opt { optVerbosity = aChange (optVerbosity opt) }++--options :: Options -> IO Options+options :: [OptDescr (Options -> IO Options)]+options = + [Option "v" ["verbose"] (NoArg (changeVerbosity (+ 1)))+ "Increases log verbosity.",+ Option "q" ["quiet"] (NoArg (changeVerbosity (flip (-) 1)))+ "Decreases log verbosity.",+ Option "" ["errors-only"] (NoArg (\opt -> return opt { eventFilter = isParseError }))+ "Show parse errors only.",+ Option "" ["unhandled"] (NoArg (\opt -> return opt { eventFilter = isUnhandled }))+ "Show parse errors only.",+ Option "p" ["printer"] (NoArg (\opt -> return opt { printing = True }))+ "Show parse errors only.",+ Option "V" ["version"] (NoArg (exitAfter showVersion ExitSuccess))+ "Print program version.",+ Option "h" ["help"] (NoArg (exitAfter showHelp ExitSuccess))+ "Prints help"+ ]++isParseError fname (PDBEvents.PDBParseError l c s) = Just $ show fname ++ ":" +++ show (PDBEvents.PDBParseError l c s)+isParseError _ _ = Nothing++isUnhandled fname (PDBEvents.PDBParseError l c s) = Just $ show fname ++ ":" ++ show (PDBEvents.PDBParseError l c s)+isUnhandled fname (PDBEvents.PDBIgnoredLine s ) = Just $ show fname ++ ":" ++ BS.unpack s+isUnhandled fname _ = Nothing+++processFile :: Options -> Bool -> String -> IO ()+processFile opts printFilename filename = (do+ performGC -- should munmap previous file and deallocate strings;+ -- no need to do it after last file, takes little before first :-).+ --input <- BS.readFile filename + input <- OpenAnyFile.readFile filename+ (Control.Monad.when printFilename $+ putStr $ concat ["File ", filename, " contains ", show (BS.length input), " characters.\n"])+ if not (printing opts)+ then parsePDBRecords filename input (processFileLoop filename opts ) ()+ else parsePDBRecords filename input (\_ e ->+ if PDBEventPrinter.isPrintable e+ then PDBEventPrinter.print System.IO.stdout e+ else return ())+ ())+ `Control.Exception.catch` + (\e -> putStrLn $ concat [show filename,+ " exception: ",+ show (e :: SomeException)])++-- Passes all events thru a filter, and prints string, if given back.+processFileLoop :: String -> Options -> () -> PDBEvents.PDBEvent -> IO ()+processFileLoop filename opts () evt = let fevt = eventFilter opts filename evt in+ case fevt of+ Just s -> putStrLn s+ Nothing -> return ()++main = do+ args <- getArgs+ let (actions, filenames, errors) = getOpt RequireOrder options args+ opts <- foldl (>>=) (return defaultOptions) actions+ let Options { optVerbosity = verbosity } = opts+ -- <- Put your code here...+ -- mapM processFile+ let show_filenames = length filenames > 1+ mapM_ (processFile opts show_filenames) filenames+ exitWith ExitSuccess++
+ tests/PrintStructureObject.hs view
@@ -0,0 +1,11 @@+{-# LANGUAGE OverloadedStrings #-}++module Main where++import System.Environment+import Bio.PDB.IO as PDBIO+import Data.ByteString.Char8 as BS++main = do [inpfname, outfname] <- getArgs+ Just structure <- PDBIO.parse $ BS.pack inpfname+ PDBIO.write structure $ BS.pack outfname