diff --git a/BioInf/MCFoldDP.hs b/BioInf/MCFoldDP.hs --- a/BioInf/MCFoldDP.hs +++ b/BioInf/MCFoldDP.hs @@ -1,6 +1,9 @@ +{-# LANGUAGE TupleSections #-} {-# LANGUAGE RecordWildCards #-} -- | +-- +-- TODO bonus system for matching "()" bracking still broken?! module BioInf.MCFoldDP where @@ -8,7 +11,7 @@ import Control.Monad import Control.Monad.ST import Data.List (find,sort) -import Data.Tuple.Select -- (sel5) +import Data.Tuple.Select (sel3) import qualified Data.Vector.Unboxed as VU import Biobase.DataSource.MCFold @@ -17,6 +20,7 @@ import Biobase.RNA.Hashes import Biobase.RNA.NucBounds import Biobase.Structure +import Biobase.Structure.Constraint import Data.PrimitiveArray import Data.PrimitiveArray.Ix @@ -24,14 +28,15 @@ -- | Folding wrapper -fold :: MotifDB -> Primary -> Tables -fold db inp = runST $ foldST db inp +fold :: MotifDB -> Primary -> Constraint -> Tables +fold db inp constr = runST $ foldST db inp constr -- | Folding in the ST monad. the number of dncm tables is the same as the -- number of known double NCMs. -foldST :: MotifDB -> Primary -> ST s Tables -foldST db inp = do +foldST :: MotifDB -> Primary -> Constraint -> ST s Tables +foldST db inp constr = do + unless (VU.length inp == (VU.length . unConstraint $ constr)) $ do error $ "foldST: inp / constr length mismatch" let n = VU.length inp -1 (sncm,sncmM) <- mkTable2 n (dncms,dncmMs) <- liftM unzip . mapM (const $ mkTable2 n) $ VU.toList knownDoubleNCM @@ -46,19 +51,18 @@ -- -- TODO check if large interior loops are ok?! Do we even want those, or -- should they come solely from NCMs? Have an option for that - let sncmIJ = fncmSingle db inp i j - -- let multiIJ = if j-i-1>0 then (VU.minimum . VU.map (\k -> mbr!(i+1,k)+mbr1!(k+1,j-1)) $ VU.enumFromN (i+1) (j-i-1)) else eInf - let multiIJ = vuminimum $ fMulti db inp i j mbr mbr1 -- if j-i-1>0 then VU.minimum . VU.map snd $ fMulti db inp i j mbr mbr1 else eInf - let interiorIJ = vuminimum $ fInterior db inp i j tmp -- VU.minimum . VU.map snd $ fInterior db inp i j tmp + let sncmIJ = fncmSingle db inp constr i j + let multiIJ = vuminimumP $ fMulti db inp constr i j mbr mbr1 + let interiorIJ = vuminimumP $ fInterior db inp constr i j tmp writeM sncmM (i,j) $ minimum [sncmIJ, multiIJ, interiorIJ] -- double NCM calculation forM (zip3 dncms dncmMs $ VU.toList knownDoubleNCM) $ \(dncm, dncmM, ((di,dj),_)) -> do let k = i+di-1 let l = j-dj+1 -- (1) Begin a stem (dNCM follows sNCM) - let dsIJ = fncmDS db inp i j k l sncm + let dsIJ = fncmDS db inp constr i j k l sncm -- (2) continue a stem (dNCM follows dNCM) - let ddIJ = VU.minimum . VU.map snd $ fncmDD db inp i j k l dncms + let ddIJ = vuminimum . VU.map snd $ fncmDD db inp constr i j k l dncms writeM dncmM (i,j) $ dsIJ `min` ddIJ -- fill helper table which makes mbr calculations ~2x faster @@ -66,9 +70,8 @@ -- fill mbr table let upIJ = if (i+1 mbr!(i,k) + tmp!(k+1,j)) $ VU.enumFromN (i+1) (j-i-1) else eInf - let mbrstemIJ = if (i+2 do let unpairedExt = extern!(i+1,j) - let stemExtExt = if (i+2 Double -> Primary -> Tables -> [(Double,Secondary)] -backtrack db delta inp (sncm,dncms,mbr,mbr1,extern) = outp where +backtrack :: MotifDB -> Double -> Primary -> Constraint -> Tables -> [(Double,D1Secondary)] +backtrack db delta inp cns@(Constraint constr) (sncm,dncms,mbr,mbr1,extern) = outp where outp - | bE >= -0.0001 = [(0,Secondary (n+1) [])] - | otherwise = map (first (bE+delta-) . second (Secondary (n+1) . sort)) $ externbt delta 0 n + | bE >= -0.0001 = [(0,f ((n+1),[]))] + | otherwise = {-filter ((<=0).fst) .-} map (first (bE-bonus+delta-) . second (f . ((n+1),) . sort)) $ externbt delta 0 n + where f :: (Int,[(Int,Int)]) -> D1Secondary + f = mkD1S externbt d i j = -- unpaired nucleotide to the left [ (e,x) @@ -127,7 +132,7 @@ | let (di,dj) = fst $ knownDoubleNCM VU.! idx , let dncm = dncms!!idx , let k = i+di-1, let l = j-dj+1 - , (nidx,bestE) <- VU.toList $ fncmDD db inp i j k l dncms + , (nidx,bestE) <- VU.toList $ fncmDD db inp cns i j k l dncms , let d' = dncm!(i,j) - bestE + d , d'>=0 , (e,x) <- dncmbt d' nidx k l @@ -137,7 +142,7 @@ | let (di,dj) = fst $ knownDoubleNCM VU.! idx , let dncm = dncms!!idx , let k = i+di-1, let l = j-dj+1 - , let bestE = fncmDS db inp i j k l sncm + , let bestE = fncmDS db inp cns i j k l sncm , let d' = dncm!(i,j) - bestE + d , d'>=0 , (e,x) <- sncmbt d' k l @@ -145,14 +150,14 @@ sncmbt d i j = [ (d',[(i,j)]) | j-i>=3 - , let bestE = fncmSingle db inp i j + , let bestE = fncmSingle db inp cns i j , let d' = sncm!(i,j) - bestE + d , d'>=0 ] ++ -- sNCM enclosing a multibranched loop [ (ey,(i,j):x++y) | j-i>3 - , (k,bestE) <- VU.toList $ fMulti db inp i j mbr mbr1 + , (k,bestE) <- VU.toList $ fMulti db inp cns i j mbr mbr1 , let d' = sncm!(i,j) - bestE + d , d'>=0 , (ex,x) <- mbrbt d' (i+1) k @@ -163,7 +168,7 @@ [ (e,(i,j):x) | idx <- [0 .. VU.length knownDoubleNCM -1] , let dncm = dncms!!idx - , ((k,l),bestE) <- VU.toList $ fInterior db inp i j dncm + , ((k,l),bestE) <- VU.toList $ fInterior db inp cns i j dncm , let d' = sncm!(i,j) - bestE + d , d'>=0 , (e,x) <- dncmbt d' idx k l @@ -221,6 +226,9 @@ ] n = VU.length inp -1 bE = extern ! (0,n) + -- bonus = bonusScore * (fromIntegral . VU.length . VU.filter (`VU.elem` bonusCC) . VU.map fst $ constr) + -- set bonus score from constraints so that if a constraint couldn't be met, we still allow the result + bonus = (fromIntegral . round $ bE / bonusScore) * bonusScore @@ -228,12 +236,15 @@ -- | singleNCM insertion -fncmSingle :: MotifDB -> Primary -> Int -> Int -> Double -fncmSingle MotifDB{..} inp i j +fncmSingle :: MotifDB -> Primary -> Constraint -> Int -> Int -> Double +fncmSingle MotifDB{..} inp cns@(Constraint constr) i j + | (fst $ constr VU.! i) == 'x' || (fst $ constr VU.! j) == 'x' = eInf | l<4 = eInf - | otherwise = maybe 0 ((!ci) . snd) $ find ((l==).fst) $ sCycles where - ci = mkHashedPrimary (minExtended,maxExtended) $ VU.slice i l inp - l = j-i+1 + | otherwise = bonus + (maybe 0 ((!ci) . snd) $ find ((l==).fst) $ sCycles) + where + ci = mkHashedPrimary (minExtended,maxExtended) $ VU.slice i l inp + l = j-i+1 + bonus = giveBonus cns i j {-# INLINE fncmSingle #-} -- | double NCM extend single NCM. We do not care that this is comparatively @@ -243,21 +254,23 @@ -- TODO otherwise case -- TODO eats another ~10% performance -fncmDS :: MotifDB -> Primary -> Int -> Int -> Int -> Int -> Table2 -> Double -fncmDS MotifDB{..} inp i j k l sncm +fncmDS :: MotifDB -> Primary -> Constraint -> Int -> Int -> Int -> Int -> Table2 -> Double +fncmDS MotifDB{..} inp cns@(Constraint constr) i j k l sncm + | (fst $ constr VU.! i) == 'x' || (fst $ constr VU.! j) == 'x' = eInf | k>=l = eInf -- in case we have one of the three known single NCMs | Just hinge <- ((di,dj),len) `lookup` dsConnect , Just ncm <- (di,dj) `lookup` dCycles - = sncm!(k,l) + hinge!(inp `VU.unsafeIndex` k, inp `VU.unsafeIndex` l) + ncm!ci + = sncm!(k,l) + hinge!(inp `VU.unsafeIndex` k, inp `VU.unsafeIndex` l) + ncm!ci + bonus -- | i==1 && j==7 && di==2 && dj==2 = error $ show (di,dj,len) -- the single NCM is larger - | otherwise = sncm!(k,l) + | otherwise = sncm!(k,l) + bonus where di = k-i+1 dj = j-l+1 len = l-k+1 ci = mkHashedPrimary (minExtended,maxExtended) $ VU.slice i di inp VU.++ VU.slice l dj inp + bonus = giveBonus cns i j {-# INLINE fncmDS #-} -- | double NCM extending another double NCM. @@ -267,18 +280,22 @@ -- TODO improve performance, eats ~66% of total time -- TODO improve: return empty vector on error, write special minimum function that has eInf on empty -fncmDD :: MotifDB -> Primary -> Int -> Int -> Int -> Int -> [Table2] -> VU.Vector (Int,Double) -fncmDD MotifDB{..} inp i j k l dncms = VU.fromList $ zipWith3 f (map fst $ VU.toList knownDoubleNCM) [0..] dncms where - f (dk,dl) idx dncm - | k+2>=l = (-1,eInf) - | Just hinge <- ((di,dj),(dk,dl)) `lookup` ddConnect - , Just ncm <- (di,dj) `lookup` dCycles - = (idx,dncm!(k,l) + hinge!(inp `VU.unsafeIndex` k, inp `VU.unsafeIndex` l) + ncm!ci) - | otherwise = (-1,eInf) - where - di = k-i+1 - dj = j-l+1 - ci = mkHashedPrimary (minExtended,maxExtended) $ VU.slice i di inp VU.++ VU.slice l dj inp +fncmDD :: MotifDB -> Primary -> Constraint -> Int -> Int -> Int -> Int -> [Table2] -> VU.Vector (Int,Double) +fncmDD MotifDB{..} inp cns@(Constraint constr) i j k l dncms + | (fst $ constr VU.! i) == 'x' || (fst $ constr VU.! j) == 'x' = VU.empty + | otherwise = VU.fromList $ zipWith3 f (map fst $ VU.toList knownDoubleNCM) [0..] dncms + where + bonus = giveBonus cns i j + f (dk,dl) idx dncm + | k+2>=l = (-1,eInf) + | Just hinge <- ((di,dj),(dk,dl)) `lookup` ddConnect + , Just ncm <- (di,dj) `lookup` dCycles + = (idx,dncm!(k,l) + hinge!(inp `VU.unsafeIndex` k, inp `VU.unsafeIndex` l) + ncm!ci + bonus) + | otherwise = (-1,eInf) + where + di = k-i+1 + dj = j-l+1 + ci = mkHashedPrimary (minExtended,maxExtended) $ VU.slice i di inp VU.++ VU.slice l dj inp {-# INLINE fncmDD #-} -- | Add one stem for "external" calculations @@ -307,9 +324,13 @@ -- -- TODO close with singleNCM -fMulti :: MotifDB -> Primary -> Int -> Int -> Table2 -> Table2 -> VU.Vector (Int,Double) -fMulti db inp i j mbr mbr1 = xs where - xs = VU.map (\k -> (k, mbr!(i+1,k) + mbr1!(k+1,j-1))) $ (VU.enumFromN (i+1) (j-i-1)) +fMulti :: MotifDB -> Primary -> Constraint -> Int -> Int -> Table2 -> Table2 -> VU.Vector (Int,Double) +fMulti db inp cns@(Constraint constr) i j mbr mbr1 + | (fst $ constr VU.! i) == 'x' || (fst $ constr VU.! j) == 'x' = VU.empty + | otherwise = xs + where + xs = VU.map (\k -> (k, mbr!(i+1,k) + mbr1!(k+1,j-1) + bonus)) $ (VU.enumFromN (i+1) (j-i-1)) + bonus = giveBonus cns i j {-# INLINE fMulti #-} -- | Connect a partial multibranched structure with a hairpin. Note that the @@ -329,33 +350,23 @@ -- -- TODO this could profit from a log-based scoring function -fInterior :: MotifDB -> Primary -> Int -> Int -> Table2 -> VU.Vector ((Int,Int),Double) -{- -fInterior MotifDB{..} inp i j dncm = VU.fromList $ ((-1,-1),eInf) : xs where - xs = [ ((k,l),dncm!(k,l)) - | k<-[i+1 .. i+10] - , l<-[j-1,j-2 .. j-10] - , k ((k,l),dncm!(k,l))) . VU.filter (\(k,l) -> k Primary -> Constraint -> Int -> Int -> Table2 -> VU.Vector ((Int,Int),Double) +fInterior MotifDB{..} inp cns@(Constraint constr) i j dncm + | (fst $ constr VU.! i) == 'x' || (fst $ constr VU.! j) == 'x' = VU.empty + | otherwise = res + where + bonus = giveBonus cns i j + res = VU.map (\(k,l) -> ((k,l),dncm!(k,l) + bonus)) . VU.filter (\(k,l) -> k Double -getScore ts = score where - t = sel5 ts - (_,(n,_)) = bounds t - score = t!(0,n) +-- * Helper functions -- | The default two-dim table @@ -374,14 +385,42 @@ -- TODO the next version of the vector library should handle "VU.map snd" -- better -vuminimum xs = VU.foldl' (\a (_,x) -> min a x) eInf xs +vuminimumP xs = VU.foldl' (\a (_,x) -> min a x) eInf xs +{-# INLINE vuminimumP #-} + +vuminimum xs = VU.foldl' (\a x -> min a x) eInf xs {-# INLINE vuminimum #-} --- vuminimum = VU.foldl' min eInf . VU.map snd --- * types and newtypes +-- * types and constants type Table2 = PrimArray (Int,Int) Double type Table4 = PrimArray (Int,Int,Int,Int) Double type Tables = (Table2,[Table2],Table2,Table2,Table2) + +bonusScore :: Double +bonusScore = (-10)^5 + +-- | Give a certain 'bonusScore' for the constraints that have been fulfilled. +-- +-- TODO should we be more lenient with constraints that would increase the +-- total energy? + +giveBonus :: Constraint -> Int -> Int -> Double +giveBonus (Constraint constr) i j + | any (=='x') [c,d] = eInf + | c == '(' && d == ')' + , k==j && l==i = bonusScore * 2 + | c == '(' || d == ')' = eInf + | c == ')' || d == '(' = eInf + | c == '>' && k>i = bonusScore + | c == '>' = eInf + | d == '<' && l=0 = bonusScore + | d == '<' = eInf + | all (`VU.elem` bonusCC) [c,d] = bonusScore * 2 + | any (`VU.elem` bonusCC) [c,d] = bonusScore + | otherwise = 0 + where + (c,k) = constr VU.! i + (d,l) = constr VU.! j diff --git a/MC-Fold-DP.cabal b/MC-Fold-DP.cabal --- a/MC-Fold-DP.cabal +++ b/MC-Fold-DP.cabal @@ -1,5 +1,5 @@ name: MC-Fold-DP -version: 0.1.0.0 +version: 0.1.0.1 author: Christian Hoener zu Siederdissen, Stephan H Bernhart, Peter F Stadler, Ivo L Hofacker copyright: Christian Hoener zu Siederdissen, 2010-2011 homepage: http://www.tbi.univie.ac.at/software/mcfolddp/ @@ -48,34 +48,49 @@ . * suboptimals: return all structures within an energy band above the ground state . + * constraint folding (fill partial structures) + . Todo: . * Boltzmann likelihood calculations . * pseudoknot calculations (currently aiming for a pknotsRG-like algorithm) +Flag fastBuild + description: build using fast GHC options (developer only!) + default: False + library build-depends: base >= 4 && < 5, vector, tuple, PrimitiveArray, - Biobase == 0.3.1.0 + Biobase == 0.3.1.1 exposed-modules: BioInf.MCFoldDP - ghc-options: - -Odph - if impl(ghc > 6.13) + if flag(fastBuild) ghc-options: - -fllvm + -O0 + else + ghc-options: + -Odph + if impl(ghc > 6.13) + ghc-options: + -fllvm executable MCFoldDP build-depends: - cmdargs == 0.6.8 + cmdargs == 0.6.8, + split main-is: MCFoldDP.hs - ghc-options: - -Odph - if impl(ghc > 6.13) + if flag(fastBuild) ghc-options: - -fllvm + -O0 + else + ghc-options: + -Odph + if impl(ghc > 6.13) + ghc-options: + -fllvm diff --git a/MCFoldDP.hs b/MCFoldDP.hs --- a/MCFoldDP.hs +++ b/MCFoldDP.hs @@ -4,72 +4,112 @@ module Main where import System.Console.CmdArgs -import Control.Monad (liftM) +import Control.Monad (liftM,when) import Text.Printf import Data.List (sortBy) import Data.Ord (comparing) +import Data.List.Split (splitEvery) import Biobase.DataSource.MCFold import Biobase.DataSource.MCFold.Import import Biobase.RNA import Biobase.Structure -import Biobase.Structure.DotBracket +import Biobase.Structure.Constraint +--import Biobase.Structure.DotBracket import BioInf.MCFoldDP -data Options = Options - { database :: FilePath - , strictInput :: Bool - , noSparseDataCorrection :: Bool - , orderSuboptimals :: Bool - , band :: Double - , oneResult :: Bool - } deriving (Data,Typeable,Show) +data Options + = Fold + { database :: FilePath + , strictInput :: Bool + , noSparseDataCorrection :: Bool + , orderSuboptimals :: Bool + , band :: Double + , oneResult :: Bool + , constraint :: Maybe String + , printConstraint :: Bool + } + | Constrain + { database :: FilePath + , strictInput :: Bool + , noSparseDataCorrection :: Bool + , orderSuboptimals :: Bool + , band :: Double + , oneResult :: Bool + , printConstraint :: Bool + } + deriving (Data,Typeable,Show) -options = Options +options = Fold { database = "./MCFOLD-DB" &= typDir &= help "path to MCFOLD-DB (default: ./MCFOLD-DB)" - , strictInput = False &= help "filter out other characters than ACGU (default: convert other characters to E and handle gracefully)" + , strictInput = False &= help "Silently ignore all sequences containing characters other than ACGU (default: convert other characters to E and handle gracefully)" , noSparseDataCorrection = False &= help "disable sparse data correction (default: enabled)" , orderSuboptimals = False &= help "sort suboptimal results by score (better scores first) (default: false)" , band = 0.1 &= help "score band above the ground state for which suboptimal results are allowed (default: 0.1)" , oneResult = False &= help "Return only one of several co-optimal structures in the backtracking phase (default: false)" - } &= summary "MCFold-DP, (c) Christian Hoener zu Siederdissen et al, 2010-2011" + , constraint = Nothing &= help "Structure constraint, follows ViennaRNA notation. Cf. \"Constraint\" major mode" + , printConstraint = False &= help "prints the constraint as second line of output. This produces three or more lines: input, constraint, MC-Fold-DP suboptimals" + } &= help "Major mode: expects input on stdin, each line is assumed to be a nucleotide sequence and suboptimal secondary structure prediction is performed." &= details [ "This program performs calculations similar to those done by MC-Fold." , "Important differences are: polynomial runtime, no pseudoknot handling," , "sparse data correction, and the possibility to gracefully handle all" , "input sequences." ] +constrain = Constrain + { + } &= help "Major mode: expects the input to be in the form or two lines each, one the sequence, one the structural constraint (followed, possibly by more characters which are ignored). Using this mode, MC-Fold-DP can postprocess ViennaRNA RNAfold input." + &= details [ "This major mode expects input in two lines each from stdin." + , "The first line is interpreted as a sequence to fold in its entirety." + , "The second line is split into words. The first word has to be the same length" + , "as the first line, other words are discarded. The complete second line is" + , "printed out by giving --printConstraint. This allows for easier comparison" + ] + +ms = [options &= auto, constrain] + &= summary "MCFold-DP, (c) Christian Hoener zu Siederdissen et al, 2010-2011" + main = do - o@Options{..} <- cmdArgs options + o <- cmdArgs $ modes ms + runMain o + +runMain o@Fold{..} = do db <- parseDir database -- TODO enable sparsity correction! - cnts <- liftM ((strictFilter strictInput) . lines) $ getContents - mapM_ (doFold o db) cnts - -strictFilter :: Bool -> [String] -> [String] -strictFilter False xs = xs -strictFilter True xs = filter (all (`elem` "ACGUacgu")) xs - --- | Executes folding a single sequence. Allows + cnts <- liftM (filter (\x -> not strictInput || strict x) . lines) $ getContents + mapM_ (doFold db oneResult band orderSuboptimals printConstraint constraint) cnts -doFold :: Options -> MotifDB -> String -> IO () -doFold Options{..} db inp = do - putStrLn inp - let pri = mkPrimary inp - let ts = fold db pri - let res = (if oneResult then take 1 else id) $ backtrack db band pri ts - mapM_ (\(e,s) -> printf "%s (%7.2f)\n" (dotbracket s) e) . (if orderSuboptimals then sortBy (comparing fst) else id) $ res - if null res then error $ "XXX " ++ inp else return () +runMain o@Constrain{..} = do + db <- parseDir database + -- TODO enable sparsity correction! + cnts <- liftM (filter (\[x,y] -> not strictInput || strict x) . splitEvery 2 . lines) $ getContents + mapM_ (\[x,y] -> doFold db oneResult band orderSuboptimals printConstraint (Just y) x) cnts +strict :: String -> Bool +strict xs = (all (`elem` "ACGUacgu")) xs +-- | Executes folding a single sequence. -runtest = do - db <- parseDir "/home/choener/tmp/mcfold/MCFOLD-DB" - {- - doFold db True 8.0 "cccaaaggg" - doFold db True 10.0 "ccccccaaagggcccaaagggggg" - -} - doFold (Options undefined False False True 2.0 False) db "UGAGUUUAUCAGCUGAUUUU" +doFold :: MotifDB -> Bool -> Double -> Bool -> Bool -> Maybe String -> String -> IO () +doFold db oneResult band orderSuboptimals printConstraint cnstr inp + | Nothing <- cnstr + = do putStrLn inp + common . mkConstraint $ replicate (length inp) '.' + | Just cns <- cnstr + , [] /= cns + , let (c:_) = words cns + , length inp == length c + = do putStrLn inp + when printConstraint $ putStrLn (cns ++ " (Constraint)") + common . mkConstraint $ c + | otherwise = error $ "doFold: " ++ show (inp,cnstr) + where + common c = do + let pri = mkPrimary inp + let ts = fold db pri c + let res' = (if oneResult then take 1 else id) $ backtrack db band pri c ts + let res = if null res' then [(0,mkD1S $ replicate (length inp) '.')] else res' + mapM_ (\(e,s) -> printf "%s (%6.2f)\n" (fromD1S s :: String) e) . (if orderSuboptimals then sortBy (comparing fst) else id) $ res