diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,4 @@ +#! /usr/bin/env runhaskell + +import Distribution.Simple +main = defaultMain diff --git a/examples/CanonicalAxes.hs b/examples/CanonicalAxes.hs new file mode 100644 --- /dev/null +++ b/examples/CanonicalAxes.hs @@ -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 " + 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 diff --git a/examples/CleanPDB.hs b/examples/CleanPDB.hs new file mode 100644 --- /dev/null +++ b/examples/CleanPDB.hs @@ -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 + diff --git a/examples/PDB2Fasta.hs b/examples/PDB2Fasta.hs new file mode 100644 --- /dev/null +++ b/examples/PDB2Fasta.hs @@ -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 + diff --git a/examples/Rg.hs b/examples/Rg.hs new file mode 100644 --- /dev/null +++ b/examples/Rg.hs @@ -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 + + diff --git a/examples/ShiftToCenter.hs b/examples/ShiftToCenter.hs new file mode 100644 --- /dev/null +++ b/examples/ShiftToCenter.hs @@ -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 diff --git a/examples/SplitModels.hs b/examples/SplitModels.hs new file mode 100644 --- /dev/null +++ b/examples/SplitModels.hs @@ -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 " + , "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 diff --git a/examples/StericClashCheck.hs b/examples/StericClashCheck.hs new file mode 100644 --- /dev/null +++ b/examples/StericClashCheck.hs @@ -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 + diff --git a/examples/Viewer.hs b/examples/Viewer.hs new file mode 100644 --- /dev/null +++ b/examples/Viewer.hs @@ -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 + + diff --git a/hPDB-examples.cabal b/hPDB-examples.cabal new file mode 100644 --- /dev/null +++ b/hPDB-examples.cabal @@ -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 + diff --git a/tests/PrintEvents.hs b/tests/PrintEvents.hs new file mode 100644 --- /dev/null +++ b/tests/PrintEvents.hs @@ -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 + + diff --git a/tests/PrintStructureObject.hs b/tests/PrintStructureObject.hs new file mode 100644 --- /dev/null +++ b/tests/PrintStructureObject.hs @@ -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