diff --git a/Data/Forest/Static/Align/Affine.hs b/Data/Forest/Static/Align/Affine.hs new file mode 100644 --- /dev/null +++ b/Data/Forest/Static/Align/Affine.hs @@ -0,0 +1,394 @@ + +module Data.Forest.Static.Align.Affine where + +import Control.Monad(forM_, unless) +import Numeric.Log +import qualified Data.Text as Text +import qualified Data.Tree as T +import qualified Data.Vector as V +import qualified Data.Vector.Fusion.Stream.Monadic as SM +import qualified Data.Vector.Generic as VG +import Text.Printf +import Unsafe.Coerce + +import ADP.Fusion.Core +import Biobase.Newick +import Data.Forest.Static (TreeOrder(..),Forest) +import Data.PrimitiveArray as PA hiding (map) +import FormalLanguage.CFG +import qualified Data.Forest.Static as F +import qualified Diagrams.TwoD.ProbabilityGrid as PG + +import ADP.Fusion.Forest.Align.RL + + + +[formalLanguage| +Verbose + +Grammar: Global +N: T -- tree +N: F -- forest +N: Z -- tree for gaps +N: Q -- sibling gap mode +N: R -- parent gap mode +N: E +T: n +S: [F,F] +[F,F] -> iter <<< [T,T] [F,F] +[F,F] -> fgap <<< [T,Z] [Q,Q] +[F,F] -> fgap <<< [Z,T] [Q,Q] +[Z,T] -> indel <<< [-,n] [R,R] +[T,Z] -> delin <<< [n,-] [R,R] +[T,T] -> align <<< [n,n] [F,F] +[F,F] -> done <<< [E,E] +[R,R] -> done <<< [E,E] +[R,R] -> pgap <<< [T,T] [R,R] +[R,R] -> pgap <<< [T,Z] [R,R] +[R,R] -> pgap <<< [Z,T] [R,R] +[Q,Q] -> done <<< [E,E] +[Q,Q] -> siter <<< [T,T] [F,F] +[Q,Q] -> sgap <<< [T,Z] [Q,Q] +[Q,Q] -> sgap <<< [Z,T] [Q,Q] +[E,E] -> finalDone <<< [e,e] +// +Outside: Labolg +Source: Global +// + +Emit: Global +Emit: Labolg +|] + +makeAlgebraProduct ''SigGlobal + +resig :: Monad m => SigGlobal m a b c d -> SigLabolg m a b c d +resig (SigGlobal gdo git gsi gal gin gde gfg gpg gsg gfi gh) = SigLabolg gdo git gsi gal gin gde gfg gpg gsg gfi gh +{-# Inline resig #-} + + +score :: Monad m => Int -> Int -> Int -> Int -> SigGlobal m Int Int Info Info +score matchSc notmatchSc delinSc affinSc = SigGlobal -- match affine deletion + { gDone = \ f -> f -- (Z:.():.()) -> 0 -- traceShow "EEEEEEEEEEEEE" 0 + , gFinalDone = \ (Z:.():.()) -> 0 -- traceShow "EEEEEEEEEEEEE" 0 + , gIter = \ t f -> t+f + , gSiter = \ t f -> t+f + , gAlign = \ (Z:.c:.b) f -> tSI glb ("ALIGN",f,c,b) $ f + if label c == label b then matchSc else notmatchSc + , gIndel = \ (Z:.():.b) f -> tSI glb ("INDEL",f,b) $ f + , gDelin = \ (Z:.c:.()) f -> tSI glb ("DELIN",f,c) $ f + , gFgap = \ t f -> tSI glb ("gap",f+t,delinSc) $ t + f + delinSc --gap open + , gPgap = \ t f -> tSI glb ("gap",f+t,affinSc) $ t + f + affinSc --gap extension + , gSgap = \ t f -> tSI glb ("gap",f+t,affinSc) $ t + f + affinSc --gap extension + , gH = SM.foldl' max (-88888) + } +{-# Inline score #-} + +part :: Monad m => Log Double -> Log Double -> Log Double -> Log Double -> Log Double -> SigGlobal m (Log Double) (Log Double) Info Info +part matchSc mismatchSc delinSc affinSc temp = SigGlobal + { gDone = \ f -> f + , gIter = \ t f -> tSI glb ("TFTFTFTFTF",t,f) $ t * f + , gFinalDone = \ (Z:.():.()) -> 1 + , gSiter = \ t f -> tSI glb ("TFTFTFTFTF",t,f) $ t * f + , gAlign = \ (Z:.a:.b) f -> tSI glb ("ALIGN",f,a,b) $ f * if label a == label b then matchSc else mismatchSc + , gIndel = \ (Z:.():.b) f -> tSI glb ("INDEL",f,b) $ f + , gDelin = \ (Z:.a:.()) f -> tSI glb ("DELIN",f,a) $ f + , gFgap = \ t f -> t * f * delinSc + , gPgap = \ t f -> t * f * affinSc + , gSgap = \ t f -> t * f * affinSc + , gH = SM.foldl' (+) 0.000000 + } +{-# Inline part #-} + + + +type Pretty' = [[T.Tree (Info,Info)]] +pretty' :: Monad m => SigGlobal m [T.Tree (Info,Info)] [[T.Tree ((Info,Info))]] Info Info +pretty' = SigGlobal + { gDone = \ f -> f -- (Z:.():.()) -> [] + , gFinalDone = \ (Z:.():.()) -> [] + , gIter = \ t f -> t++f + , gSiter = \ t f -> t++f + , gAlign = \ (Z:.a:.b) f -> [T.Node (a,b) f] + , gIndel = \ (Z:.():.b) f -> [T.Node (Info "-" 0,b) f] + , gDelin = \ (Z:.a:.()) f -> [T.Node (a,Info "-" 0) f] + , gPgap = \ t f -> t ++ f + , gSgap = \ t f -> t ++ f + , gFgap = \ t f -> t ++ f + , gH = SM.toList + } +{-# Inline pretty' #-} + + + +type Trix = TreeIxR Pre V.Vector Info I +type Tbl x = TwITbl Id Unboxed (Z:.EmptyOk:.EmptyOk) (Z:.Trix:.Trix) x +type Frst = Forest Pre V.Vector Info +type TblBt x = TwITblBt Unboxed (Z:.EmptyOk:.EmptyOk) (Z:.Trix:.Trix) Int Id Id [x] +type B = T.Tree (Info,Info) + + +-- | likelihood part +-- +--NOTE for an explanation which ITbl gets @0@ or @1@ check @runInside@, +--they have the same order requirements. + +runForward :: Frst -> Frst -> Int -> Int -> Int -> Int -> Z:.Tbl Int:.Tbl Int:.Tbl Int:.Tbl Int:.Tbl Int:.Tbl Int:.Tbl Int +runForward f1 f2 matchSc notmatchSc delinSc affinSc = let + in + mutateTablesDefault $ + gGlobal (score matchSc notmatchSc delinSc affinSc) -- costs + (iTbl 0 0 (Z:.EmptyOk:.EmptyOk) (PA.fromAssocs (Z:.minIx f1:.minIx f2) (Z:.maxIx f1:.maxIx f2) (-99999) [] )) + (iTbl 0 2 (Z:.EmptyOk:.EmptyOk) (PA.fromAssocs (Z:.minIx f1:.minIx f2) (Z:.maxIx f1:.maxIx f2) (-99998) [] )) + (iTbl 0 2 (Z:.EmptyOk:.EmptyOk) (PA.fromAssocs (Z:.minIx f1:.minIx f2) (Z:.maxIx f1:.maxIx f2) (-99997) [] )) + (iTbl 0 2 (Z:.EmptyOk:.EmptyOk) (PA.fromAssocs (Z:.minIx f1:.minIx f2) (Z:.maxIx f1:.maxIx f2) (-99996) [] )) + (iTbl 0 1 (Z:.EmptyOk:.EmptyOk) (PA.fromAssocs (Z:.minIx f1:.minIx f2) (Z:.maxIx f1:.maxIx f2) (-99995) [] )) + (iTbl 0 1 (Z:.EmptyOk:.EmptyOk) (PA.fromAssocs (Z:.minIx f1:.minIx f2) (Z:.maxIx f1:.maxIx f2) (-99994) [] )) + (iTbl 0 1 (Z:.EmptyOk:.EmptyOk) (PA.fromAssocs (Z:.minIx f1:.minIx f2) (Z:.maxIx f1:.maxIx f2) (-99993) [] )) + (node NTany $ F.label f1) + (node NTany $ F.label f2) +{-# NoInline runForward #-} + +iTbl = ITbl +{-# NoInline iTbl #-} + +-- |inside part +-- +-- NOTE : Each @ITbl@ has a big-order (1st index) and a little-order index. +-- We need these indices because we operate on unboxed tables internally +-- for efficiency reasons. +-- +-- The big-order index is for full table calculations. All @ITbl@s with +-- smaller "big order" are fully calculated before any @ITbl@s with larger +-- big orders. +-- +-- The little-order influences per-cell calculations. Cells with the same +-- index are sorted by order and those with a smaller index calculated +-- first. +-- +-- For this calculation we have the following: +-- +-- @EE@ has only one rule: @EE -> ε@, the full table is filled before any +-- other table and has big order @0@. +-- +-- Now, we only have tables with big order @1@ left. +-- +-- @TT@, @TZ@, and @ZT@ have only rules that each have at least one +-- non-empty terminal in play. @TZ@ for example has @TZ -> [n,-] RR@. +-- They are executed first for a given index, because they are guaranteed +-- to get "smaller" during recursion. Hence little order @0@. +-- +-- Finally, @FF@, @QQ@, @RR@ all combine *only* syntactic symbols on the +-- right-hand side. In addition these symbols can be empty. They need to +-- come last and have little order @1@. + +runInside :: Frst -> Frst -> Log Double -> Log Double -> Log Double -> Log Double -> Log Double -> Z:.Tbl (Log Double):.Tbl (Log Double):.Tbl (Log Double):.Tbl (Log Double):.Tbl (Log Double):.Tbl (Log Double):.Tbl (Log Double) +runInside f1 f2 matchSc notmatchSc delinSc affinSc temperature = let + in + mutateTablesDefault $ + gGlobal (part matchSc notmatchSc delinSc affinSc temperature) -- costs + (ITbl 0 0 (Z:.EmptyOk:.EmptyOk) (PA.fromAssocs (Z:.minIx f1:.minIx f2) (Z:.maxIx f1:.maxIx f2) (0.00000001) [] )) -- EE + (ITbl 1 1 (Z:.EmptyOk:.EmptyOk) (PA.fromAssocs (Z:.minIx f1:.minIx f2) (Z:.maxIx f1:.maxIx f2) (0.00000002) [] )) -- FF + (ITbl 1 1 (Z:.EmptyOk:.EmptyOk) (PA.fromAssocs (Z:.minIx f1:.minIx f2) (Z:.maxIx f1:.maxIx f2) (0.00000003) [] )) -- QQ + (ITbl 1 1 (Z:.EmptyOk:.EmptyOk) (PA.fromAssocs (Z:.minIx f1:.minIx f2) (Z:.maxIx f1:.maxIx f2) (0.00000004) [] )) -- RR + (ITbl 1 0 (Z:.EmptyOk:.EmptyOk) (PA.fromAssocs (Z:.minIx f1:.minIx f2) (Z:.maxIx f1:.maxIx f2) (0.00000005) [] )) -- TT + (ITbl 1 0 (Z:.EmptyOk:.EmptyOk) (PA.fromAssocs (Z:.minIx f1:.minIx f2) (Z:.maxIx f1:.maxIx f2) (0.00000006) [] )) -- TZ + (ITbl 1 0 (Z:.EmptyOk:.EmptyOk) (PA.fromAssocs (Z:.minIx f1:.minIx f2) (Z:.maxIx f1:.maxIx f2) (0.00000007) [] )) -- ZT + (node NTany $ F.label f1) + (node NTany $ F.label f2) +{-# NoInline runInside #-} + + + +-- outside part +type Trox = TreeIxR Pre V.Vector Info O +type OTbl x = TwITbl Id Unboxed (Z:.EmptyOk:.EmptyOk) (Z:.Trox:.Trox) x + +-- | Actual outside calculations. +-- +-- NOTE: The big and little order indices are reversed compared to +-- @runInside@. We need to evaluate the tables outside in now. + +runOutside :: Frst -> Frst -> Log Double -> Log Double -> Log Double -> Log Double -> Log Double -> Z:.Tbl (Log Double):.Tbl (Log Double):.Tbl (Log Double):.Tbl (Log Double):.Tbl (Log Double):.Tbl (Log Double):.Tbl (Log Double) -> Z:.OTbl (Log Double):.OTbl (Log Double):.OTbl (Log Double):.OTbl (Log Double):.OTbl (Log Double):.OTbl (Log Double):.OTbl (Log Double) +runOutside f1 f2 matchSc mismatchSc indelSc affinSc temperature (Z:.iE:.iF:.iQ:.iR:.iT:.iS:.iZ) + = mutateTablesDefault $ + gLabolg (resig (part matchSc mismatchSc indelSc affinSc temperature)) + (ITbl 1 1 (Z:.EmptyOk:.EmptyOk) (PA.fromAssocs (Z:.minIx f1:.minIx f2) (Z:.maxIx f1:.maxIx f2) (0.00000001) [] )) + (ITbl 0 0 (Z:.EmptyOk:.EmptyOk) (PA.fromAssocs (Z:.minIx f1:.minIx f2) (Z:.maxIx f1:.maxIx f2) (0.00000002) [] )) + (ITbl 0 0 (Z:.EmptyOk:.EmptyOk) (PA.fromAssocs (Z:.minIx f1:.minIx f2) (Z:.maxIx f1:.maxIx f2) (0.00000003) [] )) + (ITbl 0 0 (Z:.EmptyOk:.EmptyOk) (PA.fromAssocs (Z:.minIx f1:.minIx f2) (Z:.maxIx f1:.maxIx f2) (0.00000004) [] )) + (ITbl 0 1 (Z:.EmptyOk:.EmptyOk) (PA.fromAssocs (Z:.minIx f1:.minIx f2) (Z:.maxIx f1:.maxIx f2) (0.00000005) [] )) + (ITbl 0 1 (Z:.EmptyOk:.EmptyOk) (PA.fromAssocs (Z:.minIx f1:.minIx f2) (Z:.maxIx f1:.maxIx f2) (0.00000006) [] )) + (ITbl 0 1 (Z:.EmptyOk:.EmptyOk) (PA.fromAssocs (Z:.minIx f1:.minIx f2) (Z:.maxIx f1:.maxIx f2) (0.00000007) [] )) + iF + iQ + iR + iT + iS + iZ + (node NTany $ F.label f1) + (node NTany $ F.label f2) +{-# NoInline runOutside #-} + + + +-- inside part +run :: Frst -> Frst -> Int -> Int -> Int -> Int -> (Z:.Tbl Int:.Tbl Int:.Tbl Int:.Tbl Int:.Tbl Int:.Tbl Int:.Tbl Int,Int,Pretty') +run f1 f2 matchSc notmatchSc delinSc affinSc = (fwd,unId $ axiom a2, unId $ axiom b2) + where fwd@(Z:.a1:.a2:.a3:.a4:.a5:.a6:.a7) = runForward f1 f2 matchSc notmatchSc delinSc affinSc + Z:.b1:.b2:.b3:.b4:.b5:.b6:.b7 + = gGlobal ((score matchSc notmatchSc delinSc affinSc) <|| pretty') + (toBacktrack a1 (undefined :: Id a -> Id a)) + (toBacktrack a2 (undefined :: Id a -> Id a)) + (toBacktrack a3 (undefined :: Id a -> Id a)) + (toBacktrack a4 (undefined :: Id a -> Id a)) + (toBacktrack a5 (undefined :: Id a -> Id a)) + (toBacktrack a6 (undefined :: Id a -> Id a)) + (toBacktrack a7 (undefined :: Id a -> Id a)) + (node NTany $ F.label f1) (node NTany $ F.label f2) + :: Z:.TblBt B:.TblBt B:.TblBt B:.TblBt B:.TblBt B:.TblBt B:.TblBt B +{-# NoInline run #-} + + + +-- outside part + +runIO f1 f2 matchSc mismatchSc indelSc affinSc temperature = (fwd,out,unId $ axiom f) + where fwd@(Z:.e:.f:.q:.r:.t:.s:.z) = runInside f1 f2 matchSc mismatchSc indelSc affinSc temperature + out@(Z:.oet:.oft:.oqt:.ort:.ott:.ost:.ozt) = runOutside f1 f2 matchSc mismatchSc indelSc affinSc temperature fwd +{-# NoInline runIO #-} + +-- a a +-- / \ / \ +-- e d b f +-- / \ / \ +-- b c c d +-- +-- (a,a) 100 +-- / \ +-- (e,-) (-,f) (-3) (-5) +-- / \ / \ +-- (b,b) (c,-) (-,c) (d,d) 100 (-3) (-5) 100 +-- +-- +-- +-- (a,a) 100 +-- / \ +-- (e,-) (d,-) (-3) (-3) +-- / \ +-- (b,b) (-,f) 100 (-5) +-- / \ +-- (c,c) (-,d) 100 (-5) + +-- all new from here + + +test n1 n2 = do + f1 <- readFile n1 + f2 <- readFile n2 + runAlignS f1 f2 10 (-30) (-10) (-1) +{-# NoInline test #-} + + + +-- | + +runAlignScoreTrees t1 t2 matchSc notmatchSc delinSc affinSc = do + let (fwd,sc,bt') = run t1 t2 matchSc notmatchSc delinSc affinSc + let (Z:.TW (ITbl _ _ _ iet) _ :.TW (ITbl _ _ _ ift) _ :.TW (ITbl _ _ _ iqt) _ :.TW (ITbl _ _ _ irt) _ :.TW (ITbl _ _ _ itt) _ :.TW (ITbl _ _ _ ist) _ :.TW (ITbl _ _ _ izt) _) = fwd + let bt = take 1 bt' -- nub bt' + printf "Score: %10d\n" sc + forM_ bt $ \b -> do + putStrLn "" + forM_ b $ \x -> putStrLn $ T.drawTree $ fmap show x +{-# NoInline runAlignScoreTrees #-} + + + +-- | + +runAlignScoreTreesIO + ∷ PG.FillWeight + → PG.RenderChoice + → FilePath + → Frst + → Frst + → Log Double + → Log Double + → Log Double + → Log Double + → IO () +runAlignScoreTreesIO fw probFileTy probFile t1 t2 matchSc mismatchSc indelSc affinSc = do + let temperature = 1.0 + let (inn,out,_) = runIO t1 t2 matchSc mismatchSc indelSc affinSc temperature -- (t2 {F.lsib = VG.fromList [-1,-1], F.rsib = VG.fromList [-1,-1]}) + let (Z:.TW (ITbl _ _ _ iet) _ :.TW (ITbl _ _ _ ift) _ :.TW (ITbl _ _ _ iqt) _ :.TW (ITbl _ _ _ irt) _ :.TW (ITbl _ _ _ itt) _ :.TW (ITbl _ _ _ ist) _ :.TW (ITbl _ _ _ izt) _) = inn + let (Z:.TW (ITbl _ _ _ iet) _ :.TW (ITbl _ _ _ oft) _ :.TW (ITbl _ _ _ oqt) _ :.TW (ITbl _ _ _ ort) _ :.TW (ITbl _ _ _ ott) _ :.TW (ITbl _ _ _ ost) _ :.TW (ITbl _ _ _ ozt) _) = out + let (Z:.(TreeIxR frst1 lb1 _):.(TreeIxR frst2 lb2 _), Z:.(TreeIxR _ ub1 _):.(TreeIxR _ ub2 _)) = bounds oft + let ix = (Z:.TreeIxR frst1 lb1 F:.TreeIxR frst2 lb2 F) + let sc = ift ! ix + print sc + let ps = map (\(k,k1,k2) -> + let k' = unsafeCoerce k + in ( k1 + , k2 + , ((itt!k) * (ott!k') / sc) + , (maybe "-" label $ F.label t1 VG.!? k1) + , (maybe "-" label $ F.label t2 VG.!? k2) + )) [ (Z:.TreeIxR frst1 k1 T:.TreeIxR frst2 k2 T,k1,k2) | k1 <- [lb1 .. ub1 - 1], k2 <- [lb2 .. ub2 - 1] ] + -- + let gsc = map (\(k1,k2,sc,l1,l2) -> sc) ps + let fillText [] = " " + fillText xs = xs + mapM_ print gsc + let gl1 = map (\k1 -> fillText . Text.unpack $ (maybe "-" label $ F.label t1 VG.!? k1)) [lb1 .. ub1 - 1] + let gl2 = map (\k2 -> fillText . Text.unpack $ (maybe "-" label $ F.label t2 VG.!? k2)) [lb2 .. ub2 - 1] + case probFileTy of + PG.SVG -> PG.svgGridFile probFile fw PG.FSfull ub1 ub2 gl1 gl2 gsc + PG.EPS -> PG.epsGridFile probFile fw PG.FSfull ub1 ub2 gl1 gl2 gsc +{-# NoInline runAlignScoreTreesIO #-} + + + +runAlignS t1' t2' matchSc notmatchSc delinSc affinSc = do + let f x = either error (F.forestPre . map getNewickTree) $ newicksFromText x + t1 = f $ Text.pack t1' + t2 = f $ Text.pack t2' + let (fwd,sc,bt') = run t1 t2 matchSc notmatchSc delinSc affinSc + let (Z:.TW (ITbl _ _ _ iet) _ :.TW (ITbl _ _ _ ift) _ :.TW (ITbl _ _ _ iqt) _ :.TW (ITbl _ _ _ irt) _ :.TW (ITbl _ _ _ itt) _ :.TW (ITbl _ _ _ ist) _ :.TW (ITbl _ _ _ izt) _) = fwd + let bt = take 1 bt' -- nub bt' + printf "Score: %10d\n" sc +-- forM_ bt $ \b -> do +-- putStrLn "" +-- forM_ b $ \x -> putStrLn $ T.drawTree $ fmap show x +{-# NoInline runAlignS #-} + +runAlignIO fw probFileTy probFile t1' t2' matchSc mismatchSc indelSc affinSc temperature = do + let f x = either error (F.forestPre . map getNewickTree) $ newicksFromText x + t1 = f $ Text.pack t1' + t2 = f $ Text.pack t2' + let (inn,out,_) = runIO t1 t2 matchSc mismatchSc indelSc affinSc temperature -- (t2 {F.lsib = VG.fromList [-1,-1], F.rsib = VG.fromList [-1,-1]}) + let (Z:.TW (ITbl _ _ _ iet) _ :.TW (ITbl _ _ _ ift) _ :.TW (ITbl _ _ _ iqt) _ :.TW (ITbl _ _ _ irt) _ :.TW (ITbl _ _ _ itt) _ :.TW (ITbl _ _ _ ist) _ :.TW (ITbl _ _ _ izt) _) = inn + let (Z:.TW (ITbl _ _ _ iet) _ :.TW (ITbl _ _ _ oft) _ :.TW (ITbl _ _ _ oqt) _ :.TW (ITbl _ _ _ ort) _ :.TW (ITbl _ _ _ ott) _ :.TW (ITbl _ _ _ ost) _ :.TW (ITbl _ _ _ ozt) _) = out + let (Z:.(TreeIxR frst1 lb1 _):.(TreeIxR frst2 lb2 _), Z:.(TreeIxR _ ub1 _):.(TreeIxR _ ub2 _)) = bounds oft + let ix = (Z:.TreeIxR frst1 lb1 F:.TreeIxR frst2 lb2 F) + let sc = ift ! ix + print sc + let ps = map (\(k,k1,k2) -> + let k' = unsafeCoerce k + in ( k1 + , k2 + , ((itt!k) * (ott!k') / sc) + , (maybe "-" label $ F.label t1 VG.!? k1) + , (maybe "-" label $ F.label t2 VG.!? k2) + )) [ (Z:.TreeIxR frst1 k1 T:.TreeIxR frst2 k2 T,k1,k2) | k1 <- [lb1 .. ub1 - 1], k2 <- [lb2 .. ub2 - 1] ] + -- + let gsc = map (\(k1,k2,sc,l1,l2) -> sc) ps + let fillText [] = " " + fillText xs = xs + mapM_ print gsc + let gl1 = map (\k1 -> fillText . Text.unpack $ (maybe "-" label $ F.label t1 VG.!? k1)) [lb1 .. ub1 - 1] + let gl2 = map (\k2 -> fillText . Text.unpack $ (maybe "-" label $ F.label t2 VG.!? k2)) [lb2 .. ub2 - 1] + case probFileTy of + PG.SVG -> PG.svgGridFile probFile fw PG.FSfull ub1 ub2 gl1 gl2 gsc + PG.EPS -> PG.epsGridFile probFile fw PG.FSfull ub1 ub2 gl1 gl2 gsc +{-# NoInline runAlignIO #-} + diff --git a/Data/Forest/Static/Align/Linear.hs b/Data/Forest/Static/Align/Linear.hs new file mode 100644 --- /dev/null +++ b/Data/Forest/Static/Align/Linear.hs @@ -0,0 +1,322 @@ + +-- | A simple linear alignment module for forest-like structures. Currently all +-- comparisons are on Newick-like forests. + +module Data.Forest.Static.Align.Linear where + +import Control.Monad(forM_, unless) +import Numeric.Log +import qualified Data.Text as Text +import qualified Data.Tree as T +import qualified Data.Vector as V +import qualified Data.Vector.Fusion.Stream.Monadic as SM +import qualified Data.Vector.Generic as VG +import Text.Printf +import Unsafe.Coerce + +import ADP.Fusion.Core +import Biobase.Newick +import Data.Forest.Static (TreeOrder(..),Forest) +import Data.PrimitiveArray as PA hiding (map) +import FormalLanguage.CFG +import qualified Data.Forest.Static as F +import qualified Diagrams.TwoD.ProbabilityGrid as PG + +import ADP.Fusion.Forest.Align.RL + + + +[formalLanguage| +Verbose + +Grammar: Global +N: T +N: F +N: M +T: n +S: [F,F] +[F,F] -> iter <<< [T,T] [F,F] +[F,F] -> iter <<< [M,M] [F,F] +[T,T] -> indel <<< [-,n] [F,F] +[T,T] -> delin <<< [n,-] [F,F] +[M,M] -> align <<< [n,n] [F,F] +[F,F] -> done <<< [e,e] +// +Outside: Labolg +Source: Global +// + +Emit: Global +Emit: Labolg +|] + +makeAlgebraProduct ''SigGlobal + +resig :: Monad m => SigGlobal m a b c d -> SigLabolg m a b c d +resig (SigGlobal gdo git gal gin gde gh) = SigLabolg gdo git gal gin gde gh +{-# Inline resig #-} + +score :: Monad m => Int -> Int -> Int -> SigGlobal m Int Int Info Info +score matchSc notmatchSc delinSc = SigGlobal + { gDone = \ (Z:.():.()) -> 0 + , gIter = \ t f -> t+f + , gAlign = \ (Z:.a:.b) f -> f + if label a == label b then matchSc else notmatchSc + , gIndel = \ (Z:.():.b) f -> f + delinSc + , gDelin = \ (Z:.a:.()) f -> f + delinSc + , gH = SM.foldl' max (-88888) + } +{-# Inline score #-} + +part :: Monad m => Log Double -> Log Double -> Log Double -> Log Double -> SigGlobal m (Log Double) (Log Double) Info Info +part matchSc mismatchSc indelSc temp = SigGlobal + { gDone = \ (Z:.():.()) -> 1 + , gIter = \ t f -> t * f + , gAlign = \ (Z:.a:.b) f -> f * if label a == label b then matchSc else mismatchSc + , gIndel = \ (Z:.():.b) f -> f * indelSc + , gDelin = \ (Z:.a:.()) f -> f * indelSc + , gH = SM.foldl' (+) 0 + } +{-# Inline part #-} + +type Pretty = [[T.Tree (Info,Info)]] +pretty :: Monad m => SigGlobal m [T.Tree (Info,Info)] [[T.Tree ((Info,Info))]] Info Info +pretty = SigGlobal + { gDone = \ (Z:.():.()) -> [] + , gIter = \ t f -> t++f + , gAlign = \ (Z:.a:.b) f -> [T.Node (a,b) f] + , gIndel = \ (Z:.():.b) f -> [T.Node (Info "-" 0,b) f] + , gDelin = \ (Z:.a:.()) f -> [T.Node (a,Info "-" 0) f] + , gH = SM.toList + } +{-# Inline pretty #-} + + + +type Trix = TreeIxR Pre V.Vector Info I +type Tbl x = TwITbl Id Unboxed (Z:.EmptyOk:.EmptyOk) (Z:.Trix:.Trix) x +type Frst = Forest Pre V.Vector Info +type TblBt x = TwITblBt Unboxed (Z:.EmptyOk:.EmptyOk) (Z:.Trix:.Trix) Int Id Id [x] +type B = T.Tree (Info,Info) + +runForward :: Frst -> Frst -> Int -> Int -> Int -> Z:.Tbl Int :.Tbl Int:.Tbl Int +runForward f1 f2 matchSc notmatchSc delinSc = mutateTablesDefault $ + gGlobal (score matchSc notmatchSc delinSc) + (ITbl 0 1 (Z:.EmptyOk:.EmptyOk) (PA.fromAssocs (Z:.minIx f1:.minIx f2) (Z:.maxIx f1:.maxIx f2) (-99999) [] )) + (ITbl 0 0 (Z:.EmptyOk:.EmptyOk) (PA.fromAssocs (Z:.minIx f1:.minIx f2) (Z:.maxIx f1:.maxIx f2) (-99998) [] )) + (ITbl 0 0 (Z:.EmptyOk:.EmptyOk) (PA.fromAssocs (Z:.minIx f1:.minIx f2) (Z:.maxIx f1:.maxIx f2) (-99997) [] )) + (node NTany $ F.label f1) + (node NTany $ F.label f2) +{-# NoInline runForward #-} + +runInside :: Frst -> Frst -> Log Double -> Log Double -> Log Double -> Log Double -> Z:.Tbl (Log Double):.Tbl (Log Double):.Tbl (Log Double) +runInside f1 f2 matchSc mismatchSc indelSc temperature = mutateTablesDefault $ + gGlobal (part matchSc mismatchSc indelSc temperature) + (ITbl 0 1 (Z:.EmptyOk:.EmptyOk) (PA.fromAssocs (Z:.minIx f1:.minIx f2) (Z:.maxIx f1:.maxIx f2) (0.00001) [] )) + (ITbl 0 0 (Z:.EmptyOk:.EmptyOk) (PA.fromAssocs (Z:.minIx f1:.minIx f2) (Z:.maxIx f1:.maxIx f2) (0.00002) [] )) + (ITbl 0 0 (Z:.EmptyOk:.EmptyOk) (PA.fromAssocs (Z:.minIx f1:.minIx f2) (Z:.maxIx f1:.maxIx f2) (0.00003) [] )) + (node NTany $ F.label f1) + (node NTany $ F.label f2) +{-# NoInline runInside #-} + +type Trox = TreeIxR Pre V.Vector Info O +type OTbl x = TwITbl Id Unboxed (Z:.EmptyOk:.EmptyOk) (Z:.Trox:.Trox) x + +runOutside :: Frst -> Frst -> Log Double -> Log Double -> Log Double -> Log Double -> Z:.Tbl (Log Double):.Tbl (Log Double):.Tbl (Log Double) -> Z:.OTbl (Log Double):.OTbl (Log Double):.OTbl (Log Double) +runOutside f1 f2 matchSc mismatchSc indelSc temperature (Z:.iF:.iM:.iT) + = mutateTablesDefault $ + gLabolg (resig (part matchSc mismatchSc indelSc temperature)) + (ITbl 0 0 (Z:.EmptyOk:.EmptyOk) (PA.fromAssocs (Z:.minIx f1:.minIx f2) (Z:.maxIx f1:.maxIx f2) (0.00001) [] )) + (ITbl 0 1 (Z:.EmptyOk:.EmptyOk) (PA.fromAssocs (Z:.minIx f1:.minIx f2) (Z:.maxIx f1:.maxIx f2) (0.00002) [] )) + (ITbl 0 1 (Z:.EmptyOk:.EmptyOk) (PA.fromAssocs (Z:.minIx f1:.minIx f2) (Z:.maxIx f1:.maxIx f2) (0.00003) [] )) + iF + iM + iT + (node NTany $ F.label f1) + (node NTany $ F.label f2) +{-# NoInline runOutside #-} + +runS :: Frst -> Frst -> Int -> Int -> Int -> (Z:.Tbl Int :.Tbl Int:.Tbl Int, Int ,[[T.Tree (Info, Info)]] ) +runS f1 f2 matchSc notmatchSc delinSc = (fwd,unId $ axiom f, unId $ axiom fb) + where fwd@(Z:.f:.m:.t) = runForward f1 f2 matchSc notmatchSc delinSc + Z:.fb:.fm:.tb = gGlobal ((score matchSc notmatchSc delinSc) <|| pretty) (toBacktrack f (undefined :: Id a -> Id a)) (toBacktrack m (undefined :: Id a -> Id a)) (toBacktrack t (undefined :: Id a -> Id a)) + (node NTany $ F.label f1) (node NTany $ F.label f2) + :: Z:.TblBt B:.TblBt B:.TblBt B +{-# NoInline runS #-} + +runIO f1 f2 matchSc mismatchSc indelSc temperature = (fwd,out,unId $ axiom f) + where fwd@(Z:.f:.m:.t) = runInside f1 f2 matchSc mismatchSc indelSc temperature + out@(Z:.oft:.omt:.ott) = runOutside f1 f2 matchSc mismatchSc indelSc temperature fwd +{-# NoInline runIO #-} + +-- a a +-- / \ / \ +-- e d b f +-- / \ / \ +-- b c c d +-- +-- (a,a) 100 +-- / \ +-- (e,-) (-,f) (-3) (-5) +-- / \ / \ +-- (b,b) (c,-) (-,c) (d,d) 100 (-3) (-5) 100 +-- +-- +-- +-- (a,a) 100 +-- / \ +-- (e,-) (d,-) (-3) (-3) +-- / \ +-- (b,b) (-,f) 100 (-5) +-- / \ +-- (c,c) (-,d) 100 (-5) + +t11 = "a;" +t12 = "a;" +t21 = "(b,c)a;" +t22 = "(b,c)a;" +t31 = "((d,e,f)b,(z)c)a;" -- +t32 = "(((d,e)y,f)b,(c,(x)i)g)a;" -- +t41 = "d;(b)e;" -- (b,c)e;" -- '-3' +t42 = "(d)f;b;" -- b;" +t51 = "(b:1,c:1)a:1;" +t52 = "b:2;c:2;" +t61 = "((b,c)e,d)a;" +t62 = "(b,(c,d)f)a;" +t71 = "(b)a;" +t72 = "(b)a;" + + +test f1 f2 = do + n1 ← readFile f1 + n2 ← readFile f2 + runAlignS n1 n2 1 (-1) (-1) +{-# NoInline test #-} + + + +-- | This is the main entry for "main-like" calls that already provide trees. A +-- match is defined as equality on the labels. +-- +-- TODO this should be turned into the more ususual streaming-compatible +-- behaviour. + +runAlignScoreTrees + ∷ Frst + -- ^ Forest @A@ to compare + → Frst + -- ^ Forest @B@ to compare + → Int + -- ^ score for matching nodes between @A@ and @B@. + → Int + -- ^ if nodes between @A/B@ are matched up but have unequal labels, this + -- score is chosen. + → Int + -- ^ Score for an @in/del@ operation. + → IO () +runAlignScoreTrees t1 t2 matchSc notmatchSc delinSc = do + let (fwd,sc,bt') = runS t1 t2 matchSc notmatchSc delinSc + let (Z:.TW (ITbl _ _ _ ift) _ :. TW (ITbl _ _ _ imt) _ :. TW (ITbl _ _ _ itt) _) = fwd + let bt = take 1 bt' -- TODO make nice !!! nub bt' + printf "Score: %10d\n" sc + forM_ bt $ \b -> do + putStrLn "" + forM_ b $ \x -> putStrLn $ T.drawTree $ fmap show x +{-# NoInline runAlignScoreTrees #-} + + + +runAlignScoreTreesIO + ∷ PG.FillWeight + → PG.RenderChoice + → FilePath + → Frst + → Frst + → Log Double + → Log Double + → Log Double + → IO () +runAlignScoreTreesIO fw probFileTy probFile t1 t2 matchSc mismatchSc indelSc = do + let temperature = 1.0 + let (inn,out,_) = runIO t1 t2 matchSc mismatchSc indelSc temperature + let (Z:.TW (ITbl _ _ _ ift) _ :. TW (ITbl _ _ _ imt) _ :. TW (ITbl _ _ _ itt) _) = inn + let (Z:.TW (ITbl _ _ _ oft) _ :. TW (ITbl _ _ _ omt) _ :. TW (ITbl _ _ _ ott) _) = out + let (Z:.(TreeIxR frst1 lb1 _):.(TreeIxR frst2 lb2 _), Z:.(TreeIxR _ ub1 _):.(TreeIxR _ ub2 _)) = bounds oft + let ix = (Z:.TreeIxR frst1 lb1 F:.TreeIxR frst2 lb2 F) + let scift = ift ! ix + print scift + let scoft = Prelude.sum [ oft ! (Z:.TreeIxR frst1 b1 F :. TreeIxR frst2 b2 F) | b1 <- [lb1 .. ub1], b2 <- [lb2 .. ub2] ] + print scoft + let scimt = Prelude.sum [ imt ! (Z:.TreeIxR frst1 b1 T :. TreeIxR frst2 b2 T) | b1 <- [lb1 .. ub1], b2 <- [lb2 .. ub2] ] + print scimt + let scomt = Prelude.sum [ omt ! (Z:.TreeIxR frst1 b1 T :. TreeIxR frst2 b2 T) | b1 <- [lb1 .. ub1], b2 <- [lb2 .. ub2] ] + print scomt + let ps = map (\(k,k1,k2) -> + let k' = unsafeCoerce k + in ( k1 + , k2 + , ((imt!k) * (omt!k') / scift) + , (maybe "-" label $ F.label t1 VG.!? k1) + , (maybe "-" label $ F.label t2 VG.!? k2) + )) [ (Z:.TreeIxR frst1 k1 T:.TreeIxR frst2 k2 T,k1,k2) | k1 <- [lb1 .. ub1 - 1], k2 <- [lb2 .. ub2 - 1] ] + -- + let gsc = map (\(k1,k2,sc,l1,l2) -> sc) ps + let fillText [] = " " + fillText xs = xs + let gl1 = map (\k1 -> fillText . Text.unpack $ (maybe "-" label $ F.label t1 VG.!? k1)) [lb1 .. ub1 - 1] + let gl2 = map (\k2 -> fillText . Text.unpack $ (maybe "-" label $ F.label t2 VG.!? k2)) [lb2 .. ub2 - 1] + case probFileTy of + PG.SVG -> PG.svgGridFile probFile fw PG.FSfull ub1 ub2 gl1 gl2 gsc + PG.EPS -> PG.epsGridFile probFile fw PG.FSfull ub1 ub2 gl1 gl2 gsc +{-# NoInline runAlignScoreTreesIO #-} + + + +runAlignS t1' t2' matchSc notmatchSc delinSc = do + let f x = either error (F.forestPre . map getNewickTree) $ newicksFromText x + t1 = f $ Text.pack t1' + t2 = f $ Text.pack t2' + let (fwd,sc,bt') = runS t1 t2 matchSc notmatchSc delinSc + let (Z:.TW (ITbl _ _ _ ift) _ :. TW (ITbl _ _ _ imt) _ :. TW (ITbl _ _ _ itt) _) = fwd + let bt = take 1 bt' -- TODO make nice !!! nub bt' + printf "Score: %10d\n" sc + forM_ bt $ \b -> do + putStrLn "" + forM_ b $ \x -> putStrLn $ T.drawTree $ fmap show x +{-# NoInline runAlignS #-} + +runAlignIO fw probFileTy probFile t1' t2' matchSc mismatchSc indelSc temperature = do + let f x = either error (F.forestPre . map getNewickTree) $ newicksFromText x + t1 = f $ Text.pack t1' + t2 = f $ Text.pack t2' + let (inn,out,_) = runIO t1 t2 matchSc mismatchSc indelSc temperature -- (t2 {F.lsib = VG.fromList [-1,-1], F.rsib = VG.fromList [-1,-1]}) + let (Z:.TW (ITbl _ _ _ ift) _ :. TW (ITbl _ _ _ imt) _ :. TW (ITbl _ _ _ itt) _) = inn + let (Z:.TW (ITbl _ _ _ oft) _ :. TW (ITbl _ _ _ omt) _ :. TW (ITbl _ _ _ ott) _) = out + let (Z:.(TreeIxR frst1 lb1 _):.(TreeIxR frst2 lb2 _), Z:.(TreeIxR _ ub1 _):.(TreeIxR _ ub2 _)) = bounds oft + let ix = (Z:.TreeIxR frst1 lb1 F:.TreeIxR frst2 lb2 F) + let scift = ift ! ix + print scift + let scoft = Prelude.sum [ oft ! (Z:.TreeIxR frst1 b1 F :. TreeIxR frst2 b2 F) | b1 <- [lb1 .. ub1], b2 <- [lb2 .. ub2] ] + print scoft + let scimt = Prelude.sum [ imt ! (Z:.TreeIxR frst1 b1 T :. TreeIxR frst2 b2 T) | b1 <- [lb1 .. ub1], b2 <- [lb2 .. ub2] ] + print scimt + let scomt = Prelude.sum [ omt ! (Z:.TreeIxR frst1 b1 T :. TreeIxR frst2 b2 T) | b1 <- [lb1 .. ub1], b2 <- [lb2 .. ub2] ] + print scomt + let ps = map (\(k,k1,k2) -> + let k' = unsafeCoerce k + in ( k1 + , k2 + , ((imt!k) * (omt!k') / scift) + , (maybe "-" label $ F.label t1 VG.!? k1) + , (maybe "-" label $ F.label t2 VG.!? k2) + )) [ (Z:.TreeIxR frst1 k1 T:.TreeIxR frst2 k2 T,k1,k2) | k1 <- [lb1 .. ub1 - 1], k2 <- [lb2 .. ub2 - 1] ] + -- + let gsc = map (\(k1,k2,sc,l1,l2) -> sc) ps + let fillText [] = " " + fillText xs = xs + let gl1 = map (\k1 -> fillText . Text.unpack $ (maybe "-" label $ F.label t1 VG.!? k1)) [lb1 .. ub1 - 1] + let gl2 = map (\k2 -> fillText . Text.unpack $ (maybe "-" label $ F.label t2 VG.!? k2)) [lb2 .. ub2 - 1] + case probFileTy of + PG.SVG -> PG.svgGridFile probFile fw PG.FSfull ub1 ub2 gl1 gl2 gsc + PG.EPS -> PG.epsGridFile probFile fw PG.FSfull ub1 ub2 gl1 gl2 gsc +{-# NoInline runAlignIO #-} + diff --git a/Forestry.cabal b/Forestry.cabal new file mode 100644 --- /dev/null +++ b/Forestry.cabal @@ -0,0 +1,183 @@ +name: Forestry +version: 0.0.0.1 +author: Christian Hoener zu Siederdissen, Sarah Berkemer, 2017 +copyright: Christian Hoener zu Siederdissen, 2017 +homepage: https://github.com/choener/Forestry +bug-reports: https://github.com/choener/Forestry/issues +maintainer: choener@bioinf.uni-leipzig.de +category: Formal Languages, Bioinformatics +license: BSD3 +license-file: LICENSE +build-type: Simple +stability: experimental +cabal-version: >= 1.10.0 +tested-with: GHC == 8.0.2 +synopsis: Comparison of trees and forests. +description: + Wikipedia on forestry: "the science and craft of creating, + managing, using, conserving, and repairing forests and + associated resources." + . + This is a modest starting point and provides just: + . + * forest alignments + * tree alignments + . + Build this library and executable with GHC-8.0.2 + . + To enable @RNAforestry@, use @cabal new-build -frnaforestry@ or + @cabal new-configure -w compiler-path -frnaforestry && cabal + new-build@. We suggest using @cabal@ with the nix-style local + build options. + + + +Extra-Source-Files: + changelog.md + README.md + + + +flag rnaforestry + description: build the RNAforestry example application + default: False + manual: True + + + +library + build-depends: base >= 4.7 && < 5.0 + , containers + , lens >= 4.0 + , log-domain >= 0.11 + , text >= 1.0 + , vector >= 0.10 + -- + , ADPfusion == 0.5.2.* + , ADPfusionForest == 0.0.0.* + , BiobaseNewick == 0.0.0.* + , BiobaseTypes == 0.1.3.* + , ForestStructures == 0.0.0.* + , FormalGrammars == 0.3.1.* + , PrimitiveArray == 0.8.0.* + , PrimitiveArray-Pretty == 0.0.0.* + exposed-modules: + Data.Forest.Static.Align.Affine + Data.Forest.Static.Align.Linear + default-language: + Haskell2010 + default-extensions: BangPatterns + , UnicodeSyntax + , TemplateHaskell + , QuasiQuotes + , TypeOperators + , DataKinds + , MultiParamTypeClasses + , TypeFamilies + , OverloadedStrings + , FlexibleContexts + ghc-options: + -- TODO switch to O2 again! + -O2 -funbox-strict-fields + + + +executable RNAforestry + if flag(rnaforestry) + buildable: + True + build-depends: base + , cmdargs >= 0.10 + , containers + , filepath + , lens + , log-domain >= 0.10 + , streaming >= 0.1 + , streaming-bytestring >= 0.1 + , streaming-utils >= 0.1 + , text + , vector + -- + , BiobaseNewick == 0.0.0.* + , ForestStructures + , PrimitiveArray + , PrimitiveArray-Pretty + , Forestry + , ViennaRNA-extras == 0.0.0.* + , BiobaseXNA == 0.10.0.* + , BiobaseTypes + else + buildable: + False + hs-source-dirs: + src + main-is: + RNAforestry.hs + default-language: + Haskell2010 + default-extensions: BangPatterns + , DataKinds + , DeriveDataTypeable + , FlexibleContexts + , GADTs + , MultiParamTypeClasses + , OverloadedStrings + , QuasiQuotes + , RecordWildCards + , TemplateHaskell + , TypeFamilies + , TypeOperators + , UnicodeSyntax + ghc-options: + -O2 + -funbox-strict-fields + + + +test-suite properties + type: + exitcode-stdio-1.0 + main-is: + properties.hs + ghc-options: + -threaded -rtsopts -with-rtsopts=-N -O2 -funbox-strict-fields + hs-source-dirs: + tests + default-language: + Haskell2010 + default-extensions: BangPatterns + , UnicodeSyntax + build-depends: base + , QuickCheck + , tasty >= 0.11 + , tasty-quickcheck >= 0.8 + , tasty-th >= 0.1 + -- + , Forestry + + + +benchmark benchmark + default-extensions: BangPatterns + , UnicodeSyntax + build-depends: base + , criterion >= 1.0.2 + -- + , Forestry + default-language: + Haskell2010 + hs-source-dirs: + tests + main-is: + benchmark.hs + type: + exitcode-stdio-1.0 + ghc-options: + -O2 + + + +source-repository head + type: git + location: git://github.com/choener/Forestry + diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright Christian Hoener zu Siederdissen 2017 + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * 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. + + * Neither the name of Christian Hoener zu Siederdissen nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS 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 COPYRIGHT +OWNER 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/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,17 @@ +[![Build Status](https://travis-ci.org/choener/Forestry.svg?branch=master)](https://travis-ci.org/choener/Forestry) + +# Forestry: Science and craft of forests + +This library provides linear and affine tree alignment algorithms. In addition, +a very simple program @RNAforestry@ is provided. This program mimics +@RNAforester@. + + + +#### Contact + +Christian Hoener zu Siederdissen +Leipzig University, Leipzig, Germany +choener@bioinf.uni-leipzig.de +http://www.bioinf.uni-leipzig.de/~choener/ + diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/changelog.md b/changelog.md new file mode 100644 --- /dev/null +++ b/changelog.md @@ -0,0 +1,5 @@ +0.0.0.1 +------- + +- initial checkin + diff --git a/src/RNAforestry.hs b/src/RNAforestry.hs new file mode 100644 --- /dev/null +++ b/src/RNAforestry.hs @@ -0,0 +1,140 @@ + +{-# Options_GHC -fno-cse #-} + +module Main where + +import Control.Lens +import Data.ByteString.Lens +import Control.Monad (void, when) +import Data.Function ((&)) +import qualified Data.Attoparsec.ByteString.Streaming as A +import qualified Data.ByteString.Streaming.Char8 as Q +import qualified Streaming as S +import qualified Streaming.Prelude as SP +import System.Console.CmdArgs +import System.Exit (exitSuccess, exitFailure) +import qualified Data.Vector as V +import Numeric.Log + +import Biobase.Newick.Types +import Biobase.Secondary.New (parseVienna, toTree, _label) +import Biobase.Types.Structure +import BioInf.ViennaRNA.RNAfold as RNAfold +import Data.Forest.Static +import qualified Diagrams.TwoD.ProbabilityGrid as PG + +import qualified Data.Forest.Static.Align.Affine as AA +import qualified Data.Forest.Static.Align.Linear as AL + + + +-- | Which options do we provide. + +data Options + = AlignLinear + { smatch ∷ Int + , smismatch ∷ Int + , sindel ∷ Int + , probabilities ∷ Bool + , fillweight ∷ PG.FillWeight + , temperature ∷ Double + } + | AlignAffine + { smatch ∷ Int + , smismatch ∷ Int + , sindelopen ∷ Int + , sindelcont ∷ Int + , probabilities ∷ Bool + , fillweight ∷ PG.FillWeight + , temperature ∷ Double + } + deriving (Show,Data,Typeable) + +oAlignLinear = AlignLinear + { smatch = 2 &= help "match score (2)" + , smismatch = (-2) &= help "mismatch score (-2)" + , sindel = (-1) &= help "in/del score (-1)" + , probabilities = False &= help "generate probability plot (it is advisable to pipe in RNAfold input with sequence identifiers)" + , fillweight = PG.FWlog &= help "" + , temperature = 0.1 + } + +oAlignAffine = AlignAffine + { smatch = 2 + , smismatch = (-2) + , sindelopen = (-1) + , sindelcont = (-1) + , probabilities = False + , fillweight = PG.FWlog + , temperature = 0.1 + } + +main ∷ IO () +main = do + o ← cmdArgs $ modes [oAlignLinear &= auto, oAlignAffine] + case o of + AlignLinear{} → runAlignLinear o + AlignAffine{} → runAlignAffine o + +-- | Read in pairs of @RNAfold@ outputs which are then fed to the linear +-- version of the tree alignment algorithm. + +runAlignLinear ∷ Options → IO () +runAlignLinear AlignLinear{..} = do + let go ([] S.:> r) = return r + go ([(i,x)] S.:> r) = return r + go ([(i,x),(j,y)] S.:> r) = do + let t1 = mkForest x + let t2 = mkForest y + AL.runAlignScoreTrees t1 t2 smatch smismatch sindel + when probabilities $ do + let fp' = x^.sequenceID.unpackedChars ++ y^.sequenceID.unpackedChars + let fp = (++ "-" ++ show i ++ "-" ++ show j ++".eps") $ if null fp' then "prob" else fp' + let t x = Exp $ fromIntegral x / temperature + AL.runAlignScoreTreesIO + fillweight PG.EPS + fp + t1 t2 + (t smatch) (t smismatch) (t sindel) + return r + go ∷ SP.Of [(Int,RNAfold)] a → IO a + rest ← Q.getContents + & A.parsed (RNAfold.pRNAfold RNAfold.ForceRNA 37) + & void + & SP.zip (SP.each [1..]) + & S.chunksOf 2 + & S.mapped SP.toList + & S.mapsM_ go + return () + +runAlignAffine ∷ Options → IO () +runAlignAffine AlignAffine{..} = do + let go ([] S.:> r) = return r + go ([(i,x)] S.:> r) = return r + go ([(i,x),(j,y)] S.:> r) = do + let t1 = mkForest x + let t2 = mkForest y + AA.runAlignScoreTrees t1 t2 smatch smismatch sindelopen sindelcont + when probabilities $ do + let fp' = x^.sequenceID.unpackedChars ++ y^.sequenceID.unpackedChars + let fp = (++ "-" ++ show i ++ "-" ++ show j ++".eps") $ if null fp' then "prob" else fp' + let t x = Exp $ fromIntegral x / temperature + AA.runAlignScoreTreesIO + fillweight PG.EPS + fp + t1 t2 + (t smatch) (t smismatch) (t sindelopen) (t sindelcont) + return r + go ∷ SP.Of [(Int,RNAfold)] a → IO a + rest ← Q.getContents + & A.parsed (RNAfold.pRNAfold RNAfold.ForceRNA 37) + & void + & SP.zip (SP.each [1..]) + & S.chunksOf 2 + & S.mapped SP.toList + & S.mapsM_ go + return () + +mkForest ∷ RNAfold → Forest 'Pre V.Vector Info +mkForest x = forestPre [toTree (const $ Just $ Info "" 0) (Info "" 0) . either (error . show) id . parseVienna $ x^.mfe.foldedStructure.rnass] + diff --git a/tests/benchmark.hs b/tests/benchmark.hs new file mode 100644 --- /dev/null +++ b/tests/benchmark.hs @@ -0,0 +1,6 @@ + +module Main where + +main ∷ IO () +main = return () + diff --git a/tests/properties.hs b/tests/properties.hs new file mode 100644 --- /dev/null +++ b/tests/properties.hs @@ -0,0 +1,6 @@ + +module Main where + +main ∷ IO () +main = return () +