MutationOrder (empty) → 0.0.0.1
raw patch · 11 files changed
+2147/−0 lines, 11 filesdep +ADPfusiondep +ADPfusionSetdep +BiobaseXNAsetup-changed
Dependencies added: ADPfusion, ADPfusionSet, BiobaseXNA, DPutils, FormalGrammars, MutationOrder, PrimitiveArray, PrimitiveArray-Pretty, QuickCheck, ShortestPathProblems, ViennaRNA-bindings, aeson, base, bimaps, bytestring, cereal, cereal-vector, cmdargs, containers, deepseq, directory, filepath, log-domain, parallel, serialize-instances, tasty, tasty-quickcheck, tasty-th, text, unordered-containers, vector, vector-strategies, zlib
Files
- BioInf/MutationOrder.hs +411/−0
- BioInf/MutationOrder/EdgeProb.hs +140/−0
- BioInf/MutationOrder/MinDist.hs +303/−0
- BioInf/MutationOrder/RNA.hs +282/−0
- LICENSE +675/−0
- MutationOrder.cabal +166/−0
- README.md +48/−0
- Setup.hs +2/−0
- changelog.md +6/−0
- src/MutationOrder.hs +100/−0
- tests/properties.hs +14/−0
+ BioInf/MutationOrder.hs view
@@ -0,0 +1,411 @@++-- | Run all steps of the HoxCluster algorithms in order.+--+-- This will produce the following:+-- +-- 1. run the minimal distance algorithm, give the minimal distance score+-- and return all co-optimal paths+--+-- 2. run the end-probability algorithm and return the probability that+-- each node is the begin/end of a chain+--+-- 3. run the edge probability algorithm and give the probability for each+-- @from :-> to@ edge+--+-- 4. with the edge probabilities, run the maximal probability path+-- algorithm, return that probability and all co-optimal paths+--+-- TODO -Pretty should yield a structure to be given to the eps or svg+-- generator. This allows more flexibility. Does diagrams offer+-- serialization?+--+-- TODO All this should be wrapped and available as a function. not just+-- providing output files.++module BioInf.MutationOrder+ ( module BioInf.MutationOrder+ , FillWeight (..)+ , FillStyle (..)+ , ScaleFunction (..)+ ) where++import qualified Data.Vector.Unboxed as VU+import Data.Tuple (swap)+import Control.Monad (unless,forM_,when)+import Data.Bits+import Data.ByteString (ByteString)+import Data.Function (on)+import Data.List (groupBy,sortBy)+import Data.Ord (comparing)+import Numeric.Log+import qualified Data.ByteString.Char8 as BS+import qualified Data.HashMap.Strict as HM+import qualified Data.Map.Strict as M+import qualified Data.Text as T+import qualified Data.Text.IO as T+import System.Directory (doesFileExist)+import System.Exit (exitFailure)+import Text.Printf+import Control.Arrow (first,second)+import System.IO (withFile,IOMode(WriteMode),hPutStrLn,Handle)+import System.Exit (exitSuccess)++import ADP.Fusion.Term.Edge.Type (From(..),To(..))+import Data.PrimitiveArray (fromEdgeBoundaryFst, EdgeBoundary(..), (:.)(..), getBoundary)+import Data.PrimitiveArray.ScoreMatrix+import qualified Data.PrimitiveArray as PA+import Diagrams.TwoD.ProbabilityGrid+import qualified Data.Bijection.HashMap as B+import qualified ShortestPath.SHP.Edge.MinDist as SHP+import Biobase.Secondary.Diagrams (d1Distance)++import BioInf.MutationOrder.EdgeProb+import BioInf.MutationOrder.MinDist+import BioInf.MutationOrder.RNA++++runMutationOrder verbose fw fs fwdScaleFunction probScaleFunction cooptCount cooptPrint lkupFile outprefix workdb temperature equalStart [ancestralFP,currentFP] = do+ -- only run if out file(s) do not exist+ dfe <- doesFileExist (outprefix ++ ".run")+ when dfe $ do+ printf "%s.run exists, ending now!\n" outprefix+ exitSuccess+ withFile (outprefix ++ ".run") WriteMode $ \oH -> do+ printf "%s.run job started!\n" outprefix+ --+ -- Initial stuff and debug information+ --+ ancestral <- stupidReader ancestralFP+ current <- stupidReader currentFP+ lkup <- case lkupFile of {Nothing -> return Nothing; Just f -> Just <$> qlines f}+ ls <- withDumpFile oH workdb ancestral current . fst $ createRNAlandscape lkup verbose ancestral current+ let mpks = sortBy (comparing snd) . B.toList $ mutationPositions ls+ let bitToNuc = M.fromList $ map (swap . first (+1)) mpks+ let nn = length mpks+ hPrintf oH "number of mutations: %d\n" $ mutationCount ls+ hPrintf oH "\n%s\n\n" $ replicate 80 '='+ --+ -- Run co-optimal lowest energy changes+ --+ let (e,bs) = runCoOptDist fwdScaleFunction ls+ let (ecount,countcount) = runCount fwdScaleFunction ls+ -- split co-optimals into "want to print" and "want to count";+ -- @countbs@ should be possible to stream+ let (printbs,countbs) = splitAt cooptPrint bs+ hPrintf oH "Best energy gain: %10.4f\n" e+ hPrintf oH "Number of co-optimal paths: %10d\n" countcount -- ((length printbs) + (length $ take (cooptCount-cooptPrint) bs))+ forM_ printbs (T.hPutStrLn oH)+ hPrintf oH "%s\n\n" $ replicate 80 '='+ --+ -- Run @First@ probability algorithm to determine the probability for+ -- each mutation to be the initial one+ --+ -- TODO this is completely wrong, because it still starts at the+ -- ancestral sequence. We would have to start at the extant sequence.+ -- Need to later think about this. But do not use any @First@ functions+ -- now!+ {-+ hPrintf oH "Chain begin probabilities:\n"+ let fps = boundaryPartFunFirst Nothing probScaleFunction ls+ forM_ mpks $ \(mp,k) -> hPrintf oH " %6d" (mp+1)+ hPrintf oH "\n"+ forM_ fps $ \(_, Exp p) -> hPrintf oH " %6.4f" (exp p)+ hPrintf oH "\n\n"+ printf "\n"+ -}+ --+ -- Run @Last@ probability algorithm to determine the probability for+ -- each mutation to be the last one+ --+ hPrintf oH "Chain end probabilities:\n"+ let fps = boundaryPartFunLast Nothing probScaleFunction ls+ forM_ mpks $ \(mp,k) -> hPrintf oH " %6d" (mp+1)+ hPrintf oH "\n"+ forM_ (bpNormalized fps) $ \(_, Exp p) -> hPrintf oH " %6.4f" (exp p)+ hPrintf oH "\n\n%s\n\n" $ replicate 80 '='+ --printf "\n"+ --+ -- Run specialized versions of the above, restricting the first mutation+ -- to the given one. Marginalized over the last probability, and rescaled+ -- we get the first probability. Completely printed out, we get the joint+ -- probability for each @i,j@ to be @first,last@ in the chain.+ --+ hPrintf oH "Restricted chain end probabilities\n"+ let rbps = map (\(mp,k) -> (mp,k,boundaryPartFunLast (Just k) probScaleFunction ls)) mpks+ {-+ forM_ rbps $ \(mp,k,bp) -> do+ hPrintf oH "%5d %5d\n" (mp+1) k+ forM_ (bpUnnormalized bp) $ \(l,Exp p) -> hPrintf oH "%7d " (bitToNuc M.! getBoundary l)+ hPrintf oH "\n"+ forM_ (bpUnnormalized bp) $ \(l,p) -> hPrintf oH "%7.2f " (exp . ln $ p / bpTotal bp)+ hPrintf oH "\n"+ hPrintf oH "\n"+ -}+ -- collect all restricted partition function scores and prepare for+ -- normalization+ let firstlastUn = M.fromList [ ((mp+1,bitToNuc M.! getBoundary l), logp)+ | (mp,k,bp) <- rbps, (l,logp) <- bpUnnormalized bp+ ]+ let firstlastZ = Numeric.Log.sum [ bpTotal bp | (_,_,bp) <- rbps ]+ let firstlastLogP = M.map (/firstlastZ) firstlastUn+ let firstlastP = M.map (exp . ln) firstlastLogP+ -- rowMarginals gives the total probability that the mutation order+ -- begins with this mutation.+ let rowMarginals = M.mapKeysWith (+) fst firstlastP+ -- colMarginals gives the total probability that the mutation order+ -- ends with this mutation.+ let colMarginals = M.mapKeysWith (+) snd firstlastP+ hPrintf oH " "+ forM_ (M.elems bitToNuc) $ \mut -> hPrintf oH "%6d " mut+ hPrintf oH " Σ\n"+ forM_ (M.elems bitToNuc) $ \frst -> do+ hPrintf oH "%4d " frst+ forM_ (M.elems bitToNuc) $ \lst -> hPrintf oH "%6.4f " (firstlastP M.! (frst,lst))+ hPrintf oH " %6.4f\n" $ rowMarginals M.! frst+ hPrintf oH "Σ "+ forM_ (M.elems colMarginals) $ hPrintf oH "%6.4f "+ hPrintf oH "\n\n"+ hPrintf oH "divergence from proper normalization: %10.8f\n" (1 - Prelude.sum firstlastP)+ hPrintf oH "row marginal sum %10.8f\n" (Prelude.sum rowMarginals)+ hPrintf oH "col marginal sum %10.8f\n" (Prelude.sum colMarginals)+ hPrintf oH "\n%s\n\n" $ replicate 80 '='+ -- debug on+ {-+ hPrintf oH "%f\n" $ ln firstlastZ+ hPrintf oH "%s " $ replicate 10 ' '+ forM_ (M.elems bitToNuc) $ \mut -> hPrintf oH "%10d " mut+ hPrintf oH "\n"+ forM_ (M.elems bitToNuc) $ \frst -> do+ hPrintf oH "%8d " frst+ forM_ (M.elems bitToNuc) $ \lst -> hPrintf oH "%10.4f " (ln $ firstlastUn M.! (frst,lst))+ hPrintf oH "\n"+ hPrintf oH "\n"+ hPrintf oH "%f\n" $ ln firstlastZ+ hPrintf oH "%s " $ replicate 10 ' '+ forM_ (M.elems bitToNuc) $ \mut -> hPrintf oH "%10d " mut+ hPrintf oH "\n"+ forM_ (M.elems bitToNuc) $ \frst -> do+ hPrintf oH "%8d " frst+ forM_ (M.elems bitToNuc) $ \lst -> hPrintf oH "%10.4f " ((ln $ firstlastUn M.! (frst,lst)) - ln firstlastZ)+ hPrintf oH "\n"+ hPrintf oH "\n"+ -}+ -- debug off+ -- debug on+ -- calculate first weight, unnormalized+ -- let firstUn = M.fromList [ ]+ -- debug off+ --+ --+ -- Run edge probability Inside/Outside calculations. These take quite+ -- a while longer.+ --+ let (ibs,eps) = edgeProbPartFun probScaleFunction ls+ hPrintf oH "pairwise next mutation probabilities:\n\n"+ hPrintf oH " "+ forM_ mpks $ \(mp,k) -> hPrintf oH " %6d" k+ hPutStrLn oH ""+ hPrintf oH " "+ forM_ mpks $ \(mp,k) -> hPrintf oH " %6d" (mp+1)+ hPutStrLn oH ""+ forM_ (zip (groupBy ((==) `on` (fromEdgeBoundaryFst . fst)) eps) mpks) $ \(rps,(mp,k)) -> do+ let (eb,_) = head rps+ hPrintf oH "%3d %3d" k (mp+1)+ forM_ rps $ \(eb,Exp p) -> hPrintf oH (" %6.4f") (exp p)+ hPrintf oH " %6.4f" (Prelude.sum $ map (exp . ln . snd) rps)+ hPrintf oH "\n"+ let colSums = M.fromListWith (+) [ (c,p) | ((_ :-> c),p) <- eps ]+ hPrintf oH " Σ "+ forM_ (M.toList colSums) $ \(c,Exp p) -> hPrintf oH (" %6.4f") (exp p)+ hPutStrLn oH "\n"+ gridFile [SVG,EPS] (outprefix ++ "-edge") fw fs nn nn (map (show . (+1) . fst) mpks) (map (show . (+1) . fst) mpks) (map snd eps)+ hPrintf oH "\n%s\n\n" $ replicate 80 '='+ --+ -- Generate the path with maximal edge probability+ --+ {-+ let eprobsFirst = edgeProbScoreMatrix ls (Prelude.map (Exp . log) $ M.elems colMarginals) eps+ let (Exp maxprob,mpbt) = SHP.runMaxEdgeProbFirst eprobsFirst+ hPrintf oH "Maximal Edge Log-Probability Sum: %6.4f with at least %d co-optimal paths\n" maxprob (length $ take cooptCount mpbt)+ hPutStrLn oH "first mutation to extant species\n"+ forM_ (take cooptPrint mpbt) $ \bt -> do+ let extractMut (SHP.BTnode (_:.To n)) = n+ extractMut (SHP.BTedge (From ff:.To tt)) = ff+ let mutationOrder = tail $ scanl (\set mut -> set `setBit` extractMut mut) zeroBits (reverse bt)+ let prettyPrint mut k = do+ let rna = rnas ls HM.! mut+ hPrintf oH " %3s %s\n %s MFE %6.4f\n %s CNT %6.4f\n"+ (maybe "anc" (show . (+1) . fst . (!!) mpks) k)+ (BS.unpack $ primarySequence rna)+ (BS.unpack $ mfeStructure rna)+ (mfeEnergy rna)+ (BS.unpack $ centroidStructure rna)+ (centroidEnergy rna)+ hPutStrLn oH $ replicate 8 ' ' ++ (take (BS.length $ primarySequence rna) . concat $ zipWith (\xs x -> xs ++ show x) (repeat $ " . ") (drop 1 $ cycle [0..9]))+ prettyPrint zeroBits Nothing+ forM_ (zip (reverse bt) mutationOrder) $ \case+ (SHP.BTnode (_:.To n),mut) -> prettyPrint mut (Just n)+ (SHP.BTedge (From ff:.To tt),mut) -> prettyPrint mut (Just ff)+ hPutStrLn oH ""+ hPutStrLn oH ""+ -}+ -- the rowMarginals hold the probabily to begin with a mutation. Since+ -- @Last@ goes from first to last mutation, this is what we need.+ let eplStartWeight = if equalStart+ then Prelude.map (const 1) $ M.elems rowMarginals+ else Prelude.map (Exp . log) $ M.elems rowMarginals+ let eprobsLast = edgeProbScoreMatrix ls eplStartWeight eps+ --print eprobsLast+ --print $ PA.assocs $ scoreMatrix eprobsLast+ let (Exp maxprobLast,lastLogProbs,mpbtLast') = SHP.runMaxEdgeProbLast eprobsLast+ let mpbtLast = map reverse mpbtLast'+ --print maxprobLast+ --print lastLogProbs+ --mapM_ print $ mpbtLast+ hPrintf oH "Fraction of optimal choice for each final mutation:\n"+ forM_ lastLogProbs $ \(PA.Boundary b, _) -> hPrintf oH " %6d" $ bitToNuc M.! b+ hPrintf oH "\n"+ forM_ lastLogProbs $ \(_, p) -> hPrintf oH " %6.4f" $ exp $ ln (p / Exp maxprobLast)+ hPrintf oH "\n\n"+ hPrintf oH "Maximal edge log-probability sum: %6.4f (P = %10.8f) with at least %d co-optimal paths\n" maxprobLast (exp maxprobLast) (length $ take cooptCount mpbtLast)+ hPutStrLn oH "(first mutation to extant species)\n"+ forM_ (take cooptPrint mpbtLast) $ \bt -> do+ let extractMut (SHP.BTnode (_:.To n)) = n+ extractMut (SHP.BTedge (From ff:.To tt)) = tt+ let mutationOrder = tail $ scanl (\set mut -> set `setBit` extractMut mut) zeroBits bt+ let prettyPrint mut k = do+ let rna = rnas ls HM.! mut+ hPrintf oH " %3s %s\n %s MFE %6.4f\n %s CNT %6.4f\n"+ (maybe "anc" (show . (+1) . fst . (!!) mpks) k)+ (BS.unpack $ primarySequence rna)+ (BS.unpack $ mfeStructure rna)+ (mfeEnergy rna)+ (BS.unpack $ centroidStructure rna)+ (centroidEnergy rna)+ hPutStrLn oH $ replicate 8 ' ' ++ (take (BS.length $ primarySequence rna) . concat $ zipWith (\xs x -> xs ++ show x) (repeat $ " . ") (drop 1 $ cycle [0..9]))+ prettyPrint zeroBits Nothing+ forM_ (zip bt mutationOrder) $ \case+ (SHP.BTnode (_:.To n),mut) -> prettyPrint mut (Just n)+ (SHP.BTedge (From ff:.To tt),mut) -> prettyPrint mut (Just tt)+ hPutStrLn oH ""+ hPutStrLn oH ""+ let meaOrder =+ let go = \case SHP.BTnode (_:.To n) -> n+ SHP.BTedge (From ff:.To tt) -> tt+ in map go $ concat $ take 1 mpbtLast+ let meaAnno = map (\k -> map (show . (+1) . fst) mpks !! k) meaOrder+ let meaEps = [ (ee !! k) !! l | let ee = groupBy ((==) `on` (fromEdgeBoundaryFst . fst)) eps, k <- meaOrder, l <- meaOrder ]+ gridFile [SVG,EPS] (outprefix ++ "-edge-meaorder") fw fs nn nn meaAnno meaAnno (map snd meaEps)+ --print eps+ --print meaEps+ {-+ let eprobsLast = edgeProbScoreMatrix ls (Prelude.map (Exp . log) $ M.elems rowMarginals) eps+ let (Exp maxprobLast,lastLogProbs,mpbtLast) = SHP.runMaxEdgeProbLast eprobsLast+ print maxprobLast+ print $ map (\(k,Exp p) -> (k,exp $ p - maxprobLast)) lastLogProbs+ mapM_ print $ concat $ take 2 mpbtLast+ let eprobsLast = edgeProbScoreMatrix ls (Prelude.map (Exp . log) $ M.elems colMarginals) eps+ let (Exp maxprobLast,lastLogProbs,mpbtLast) = SHP.runMaxEdgeProbLast eprobsLast+ print maxprobLast+ print $ map (\(k,Exp p) -> (k,exp $ p - maxprobLast)) lastLogProbs+ mapM_ print $ concat $ take 2 mpbtLast+ -}+{-# NoInline runMutationOrder #-}++posScaled :: Double -> Double -> ScaleFunction -> ScaleFunction+posScaled l s = scaleByFunction go where+ go d | d >= l = d ** s+ | otherwise = d+ {-# Inline go #-}+{-# Inlinable posScaled #-}++-- | Basepair distance++basepairDistanceMFE :: ScaleFunction+basepairDistanceMFE frna trna = fromIntegral $ d1Distance (mfeD1S frna) (mfeD1S trna)++basepairDistanceCentroid :: ScaleFunction+basepairDistanceCentroid frna trna = fromIntegral $ d1Distance (centroidD1S frna) (centroidD1S trna)++-- | Scale function for normal mfe delta energies++mfeDelta :: ScaleFunction+mfeDelta frna trna = mfeEnergy trna - mfeEnergy frna+{-# Inlinable mfeDelta #-}++-- | Scale function for normal centroid delta energies++centroidDelta :: ScaleFunction+centroidDelta frna trna = centroidEnergy trna - centroidEnergy frna+{-# Inlinable centroidDelta #-}++-- | Square positive "contributions", making bad moves more unlikely++squaredPositive :: ScaleFunction -> ScaleFunction+squaredPositive sf = scaleByFunction sp sf where+ sp d+ | d > 0 = d * d+ | otherwise = d+ {-# Inline sp #-}+{-# Inlinable squaredPositive #-}++-- | Scale by temperature (for probability stuff)++scaleTemperature :: Double -> ScaleFunction -> ScaleFunction+scaleTemperature t sf = scaleByFunction (/t) sf+{-# Inlinable scaleTemperature #-}++scaleByFunction f sf = \frna trna ->+ let d = sf frna trna+ in f d+{-# Inlinable scaleByFunction #-}++-- | Basepair distance++-- | Stupid fasta reader++stupidReader :: FilePath -> IO ByteString+stupidReader fp = do+ inp <- BS.lines <$> BS.readFile fp+ let xs = filter (\x -> not (BS.null x) && BS.head x /= '>') inp+ return $ BS.concat xs++-- | @withDumpFile@ is like @idIO :: a -> IO a@ in that it returns the data+-- we give to the function. However, in case the dump file exists, we read+-- it and return its contents, instead of recalculating. If it does not+-- exist, we dump the data in addition to returning it. This forces the+-- @Landscape@.++withDumpFile+ :: Handle+ -> FilePath+ -- ^ The path we store the serialized and compressed dump in+ -> ByteString+ -- ^ ancestral / origin sequence+ -> ByteString+ -- ^ destination sequence+ -> Landscape+ -- ^ the element which is to be serialized in the dump, or which would be+ -- the data in the dump+ -> IO Landscape+ -- ^ the data we put in, but maybe taken from the dump file+withDumpFile oH fp ancestral current l = do+ dfe <- doesFileExist fp+ if dfe then do+ hPrintf oH "using database %s to load sequence information\n" fp+ ls <- fromFileJSON fp+ -- now we check if we have a sane DB file+ unless (landscapeOrigin ls == ancestral && landscapeDestination ls == current) $ do+ hPutStrLn oH "ancestral or target sequence do not match those stored in the work database"+ hPutStrLn oH $ "given ancestral: " ++ BS.unpack ancestral+ hPutStrLn oH $ "DB ancestral: " ++ (BS.unpack $ landscapeOrigin ls)+ hPutStrLn oH $ "given current: " ++ BS.unpack current+ hPutStrLn oH $ "DB current: " ++ (BS.unpack $ landscapeDestination ls)+ exitFailure+ return ls+ else do+ hPrintf oH "database %s does not exist! Folding all intermediate structures. This may take a while!\n" fp+ toFileJSON fp l+ return l+
+ BioInf/MutationOrder/EdgeProb.hs view
@@ -0,0 +1,140 @@++module BioInf.MutationOrder.EdgeProb where++import Control.Arrow (second)+import Control.Monad (forM_)+import Data.List (nub,sort)+import Data.Text (Text,unpack)+import Data.Vector.Unboxed (Unbox)+import Numeric.Log+import qualified Data.Text as T+import qualified Data.Vector.Fusion.Stream.Monadic as SM+import Text.Printf+import qualified Data.Vector as V+import qualified Data.HashMap.Strict as HM+import Data.Bits+import Debug.Trace++import qualified Data.Bijection.HashMap as B+import ADP.Fusion.Core+import ADP.Fusion.EdgeBoundary+import ADP.Fusion.Set1+import Data.PrimitiveArray hiding (toList)+import Data.PrimitiveArray.ScoreMatrix+import Diagrams.TwoD.ProbabilityGrid+import FormalLanguage+import ShortestPath.SHP.Grammar.EdgeProbIO+import Data.Vector.Generic.Unstream++import BioInf.MutationOrder.RNA+import BioInf.MutationOrder.MinDist (ScaleFunction(..))++++-- | Before using @aInside@ the @ScoreMat@ needs to be scaled+--+-- TODO the @Edge@ needs to be an @EdgeWithActive@ to get the active bits+-- on the left in the set.++aInside :: Monad m => ScaleFunction -> Landscape -> SigEdgeProb m (Log Double) (Log Double) (Int:.From:.To) (Int:.To)+aInside scaled Landscape{..} = SigEdgeProb+ { edge = \x (fset:.From f:.To t) ->+ let frna = rnas HM.! (BitSet fset)+ trna = rnas HM.! (BitSet fset `xor` bit t)+ res' = scaled frna trna+ res = Exp . negate $ res'+ in+#ifdef ADPFUSION_DEBUGOUTPUT+ traceShow ("Edge",(BitSet fset,f,t),frna,trna,' ',res',res,x,x*res) $+#endif+ x * res+ , mpty = \() ->+#ifdef ADPFUSION_DEBUGOUTPUT+ traceShow "empty" $+#endif+ 1+ , node = \x (nset:.To n) ->+ let frna = rnas HM.! (BitSet nset)+ trna = rnas HM.! (BitSet nset `xor` bit n)+ res' = scaled frna trna+ res = Exp . negate $ res'+ in+#ifdef ADPFUSION_DEBUGOUTPUT+ traceShow ("Node",n,frna,trna,' ',res',res,x,res*x) $+#endif+ x * res+ , fini = \l (fset:.From f:.To t) r ->+ let frna = rnas HM.! (BitSet fset)+ trna = rnas HM.! (BitSet fset `xor` bit t)+ res' = scaled frna trna+ res = Exp . negate $ res'+ in+#ifdef ADPFUSION_DEBUGOUTPUT+ traceShow ("Fini",(BitSet fset,f,t),frna,trna,' ',res',res,l,r,l*r*res) $+#endif+ l * r * res+ , h = SM.foldl' (+) 0+-- , h = \s -> do v :: V.Vector (Log Double) <- streamToVectorM s+-- return $ Numeric.Log.sum v+ } where !n = fromIntegral mutationCount+{-# Inline aInside #-}++++type TF1 x = TwITbl Id Unboxed EmptyOk (BS1 Last I) x+type TL1 x = TwITbl Id Unboxed EmptyOk (BS1 Last O) x+type EB x = TwITbl Id Unboxed EmptyOk (EdgeBoundary C) x++++-- | Extract the individual partition scores.++edgeProbPartFun :: ScaleFunction -> Landscape -> ([(Boundary Last I, Log Double)], [(EdgeBoundary C, Log Double)])+edgeProbPartFun scaled landscape =+ let n = mutationCount landscape+ (Z:.sF:.sL:.sZ) = mutateTablesST $ gEdgeProb (aInside scaled landscape)+ (ITbl 0 0 EmptyOk (fromAssocs (BS1 0 (-1)) (BS1 (2^n-1) (Boundary $ n-1)) 0 []))+ (ITbl 1 0 EmptyOk (fromAssocs (BS1 0 (-1)) (BS1 (2^n-1) (Boundary $ n-1)) 0 []))+ (ITbl 2 0 EmptyOk (fromAssocs (0 :-> 0) (0 :-> (n-1)) 0 []))+ EdgeWithSet+ Singleton+ :: Z:.TF1 (Log Double):.TL1 (Log Double):.EB (Log Double)+ TW (ITbl _ _ _ pf ) _ = sZ+ TW (ITbl _ _ _ lkF) _ = sF+ TW (ITbl _ _ _ lkL) _ = sL+ bs' = assocs pf+ pssum = (Numeric.Log.sum $ Prelude.map snd bs') / (fromIntegral n - 1)+ ibs' = [ (Boundary b, p) | b <- [0..n] , let p = lkF ! (BS1 (2^n-1) (Boundary b)) ]+ obs' = [ (Boundary b, p) | b <- [0..n] , let p = lkF ! (BS1 (2^n-1) (Boundary b)) ]+ ibssum = Numeric.Log.sum $ Prelude.map snd ibs'+ obssum = Numeric.Log.sum $ Prelude.map snd ibs'+ ibs = Prelude.map (second (/ibssum)) ibs'+ bs = Prelude.map (second (/ibssum)) bs'+ in+#ifdef ADPFUSION_DEBUGOUTPUT+ traceShow (assocs lkF) $+ traceShow (assocs lkL) $+ traceShow (assocs pf) $+ traceShow (ibssum,obssum) $+ traceShow (bs',pssum,bs) $+#endif+ (ibs,bs)+{-# NoInline edgeProbPartFun #-}++-- | Turn the edge probabilities into a score matrix.++edgeProbScoreMatrix :: Landscape -> [Log Double] -> [(EdgeBoundary C, Log Double)] -> ScoreMatrix (Log Double)+edgeProbScoreMatrix Landscape{..} begs xs' = ScoreMatrix m nprobs names names+ where m = fromAssocs l h 0 xs+ l = (Z:.0:.0)+ h = (Z:.maximum [f | (f :-> _,_) <- xs']:.maximum [t | (_ :-> t,_) <- xs'])+ (Z:._:.hh) = h+ xs = [ ((Z:.f:.t),p) | (f :-> t, p) <- xs' ]+ (_,Z:._:.n) = bounds m+ names = V.fromList [ T.pack . show . (+1)+ . maybe (error "MutationOrder.EdgeProb.edgeProbScoreMatrix") id+ $ B.lookupR mutationPositions k | k <- [0..n]+ ]+ rowsums = HM.fromListWith (+) [ (f,p) | (f :-> t, p) <- xs' ]+ nprobs = fromAssocs 0 hh 1 $ zip [0..hh] begs -- [ (f,1-p) | f <- [0..hh], let p = rowsums HM.! f ]+
+ BioInf/MutationOrder/MinDist.hs view
@@ -0,0 +1,303 @@++-- | Calculate minimum-distance Hamiltonian Shortest Paths and+-- probabilities for starting nodes.+--+-- NOTE: We explicitly model starting nodes. For symmetrical distance+-- matrices, this reports begin/end probabilities. For asymmetrical+-- distance matrices, a second instances with @Last@ instead of @First@+-- boundary should be created to calculate begin/end probabilities+-- separately.++module BioInf.MutationOrder.MinDist where++import Control.Arrow (second)+import Control.Monad (forM_)+import Data.Bits+import Data.Data (Data)+import Data.List (nub,sort)+import Data.Text (Text)+import Data.Typeable(Typeable)+import Debug.Trace+import Numeric.Log+import qualified Data.ByteString.Char8 as BS+import qualified Data.HashMap.Strict as HM+import qualified Data.Text as T+import qualified Data.Vector.Fusion.Stream.Monadic as SM+import Text.Printf+import qualified Data.Vector.Unboxed as VU+import Data.Maybe (fromJust)+import qualified Data.Map.Strict as MS++import qualified Data.Bijection.HashMap as B+import ADP.Fusion.Core+import ADP.Fusion.Set1+import ADP.Fusion.Unit+import Data.PrimitiveArray hiding (toList,map)+import FormalLanguage+import ShortestPath.SHP.Grammar.MinDist++import BioInf.MutationOrder.RNA++++-- | Given the 'RNA' we come from and the 'RNA' we mutate into, derive the+-- gain or loss by a scaling function.++type ScaleFunction = RNA -> RNA -> Double++-- | Minimal distance algebra+--+-- TODO The two Ints are the indices of the nodes and could be replaced?++aMinDist :: Monad m => ScaleFunction -> Landscape -> SigMinDist m Double Double (Int:.From:.To) (Int:.To)+aMinDist scaled Landscape{..} = SigMinDist+ { edge = \x (fset:.From f:.To t) -> let frna = rnas HM.! (BitSet fset)+ trna = rnas HM.! (BitSet fset `setBit` f `setBit` t)+ in -- traceShow (BitSet fset, BitSet fset `setBit` f `setBit` t) $+ -- x + scaleFunction scaled (centroidEnergy trna - centroidEnergy frna)+ x + scaled frna trna+ , mpty = \() -> 0+ , node = \(nset:.To n) ->+ let frna = rnas HM.! (BitSet 0)+ trna = rnas HM.! (BitSet 0 `setBit` n)+ in -- scaleFunction scaled $ centroidEnergy trna - centroidEnergy frna+ scaled frna trna+ , fini = id+ , h = SM.foldl' min 999999+ }+{-# Inline aMinDist #-}++-- | Fused co-optimal counter!+--+-- TODO for now, @Int@ is assumed to be big enough...++aMinDistCount :: Monad m => ScaleFunction -> Landscape -> SigMinDist m (Double,Int) (Double,Int) (Int:.From:.To) (Int:.To)+aMinDistCount scaled Landscape{..} = SigMinDist+ { edge = \x (fset:.From f:.To t) -> let frna = rnas HM.! (BitSet fset)+ trna = rnas HM.! (BitSet fset `setBit` f `setBit` t)+ in -- traceShow (BitSet fset, BitSet fset `setBit` f `setBit` t, f,t, fst x, scaled frna trna) $+ -- x + scaleFunction scaled (centroidEnergy trna - centroidEnergy frna)+ (fst x + scaled frna trna, snd x)+ , mpty = \() -> (0,1)+ , node = \(nset:.To n) ->+ let frna = rnas HM.! (BitSet 0)+ trna = rnas HM.! (BitSet 0 `setBit` n)+ in -- scaleFunction scaled $ centroidEnergy trna - centroidEnergy frna+ (scaled frna trna,1)+ , fini = id+ , h = \xs -> do cntr <- SM.foldl' (\m (k,c) -> MS.insertWith (+) k c m) MS.empty xs+ -- traceShow cntr .+ return $ maybe (999999,0) fst $ MS.minViewWithKey cntr+ }+{-# Inline aMinDistCount #-}++-- | Sum over all states and collapse into boundary unscaled weights.++aInside :: Monad m => Maybe Int -> ScaleFunction -> Landscape -> SigMinDist m (Log Double) (Log Double) (Int:.From:.To) (Int:.To)+aInside restrictStartNode scaled Landscape{..} = SigMinDist+ { edge = \x (fset:.From f:.To t) -> let frna = rnas HM.! (BitSet fset)+ trna = rnas HM.! (BitSet fset `xor` bit t)+ res' = Exp . negate $ scaled frna trna+ res = x * res'+ in+ -- traceShow ("edge",fset,f,t,frna,trna,scaled frna trna, res', res) $+ maybe res (\k -> if k==t then 0 else res) restrictStartNode+ , mpty = \() -> 1+ , node = \(nset:.To n) ->+ let frna = rnas HM.! (BitSet 0)+ trna = rnas HM.! (BitSet 0 `xor` bit n)+ res = Exp . negate $ scaled frna trna+ in+ -- traceShow ("node",nset,n, frna, trna, scaled frna trna, res) $+ maybe res (\k -> if k==n then res else 0) restrictStartNode+ , fini = id+ , h = SM.foldl' (+) 0+ }+{-# Inline aInside #-}++-- | This should give the correct order of nodes independent of the+-- underlying @Set1 First@ or @Set1 Last@ because the @(From:.To)@ system+-- is agnostic over these.+--+-- TODO Use text builder++aPretty :: Monad m => ScaleFunction -> Landscape -> SigMinDist m Text [Text] (Int:.From:.To) (Int:.To)+aPretty scaled Landscape{..} = SigMinDist+ { edge = \x (fset:.From f:.To t) -> let frna = rnas HM.! (BitSet fset)+ trna = rnas HM.! (BitSet fset `setBit` f `setBit` t)+ eM = mfeEnergy trna - mfeEnergy frna+ eC = centroidEnergy trna - centroidEnergy frna+ eS = scaled frna trna+ f' = fromJust $ B.lookupR mutationPositions f+ t' = fromJust $ B.lookupR mutationPositions t+ in T.concat [x, showMut frna trna t' eM eC eS]+ , mpty = \() -> ""+ , node = \(nset:.To n) ->+ let+ frna = rnas HM.! (BitSet 0)+ trna = rnas HM.! (BitSet 0 `setBit` n)+ n' = fromJust $ B.lookupR mutationPositions n+ eM = mfeEnergy trna - mfeEnergy frna+ eC = centroidEnergy trna - centroidEnergy frna+ eS = scaled frna trna+ in T.concat [showHdr frna n', showMut frna trna n' eM eC eS]+ , fini = id+ , h = SM.toList+ } where+ showHdr frna n = T.concat+ [ T.pack $ printf "mutation mfe centr scfun "+ , T.pack $ VU.toList $ VU.replicate (BS.length $ primarySequence frna) ' ' VU.// (map (,'v') . sort . map fst $ B.toList mutationPositions)+ , T.pack $ "\n" ++ replicate 38 ' '+ , T.pack . take (BS.length $ primarySequence frna) . concat $ zipWith (\xs x -> xs ++ show x) (repeat $ " . ") (drop 1 $ cycle [0..9])+ , "\n"+ , T.pack $ printf "ancestral %5.1f %5.1f " (mfeEnergy frna) (centroidEnergy frna)+ , T.pack $ BS.unpack $ primarySequence frna+ , "\n"+ ]+ showMut frna trna n eM eC eS = T.concat+ [ T.pack $ printf "%5d %5.1f %5.1f %5.1f " (n+1) eM eC eS+ , T.pack . BS.unpack $ primarySequence trna+ , "\n"+ ]+{-# Inline aPretty #-}++-- | Count co-optimals++aCount :: Monad m => Landscape -> SigMinDist m Integer [Integer] (Int:.From:.To) (Int:.To)+aCount Landscape{..} = SigMinDist+ { edge = \x (fset:.From f:.To t) -> x+ , mpty = \() -> 1+ , node = \n -> 1+ , fini = id+ , h = \xs -> SM.foldl' (+) 0 xs >>= \x -> return [x]+ }+{-# Inline aCount #-}++++type TS1 x = TwITbl Id Unboxed EmptyOk (BS1 First I) x+type U x = TwITbl Id Unboxed EmptyOk (Unit I) x+type PF x = TwITbl Id Unboxed EmptyOk (Boundary First I) x++type TS1L x = TwITbl Id Unboxed EmptyOk (BS1 Last I) x+type PFL x = TwITbl Id Unboxed EmptyOk (Boundary Last I) x++type BT1 x b = TwITblBt Unboxed EmptyOk (BS1 First I) x Id Id b+type BTU x b = TwITblBt Unboxed EmptyOk (Unit I) x Id Id b++type BT1L x b = TwITblBt Unboxed EmptyOk (BS1 Last I) x Id Id b++++-- | Run the minimal distance algebra.+--+-- This produces one-boundary sets. Meaning that for each boundary we get+-- the total distance within the set.++forwardMinDist1 :: ScaleFunction -> Landscape -> Z:.TS1L Double:.U Double+forwardMinDist1 scaleFunction landscape =+ let n = mutationCount landscape+ in mutateTablesST $ gMinDist (aMinDist scaleFunction landscape)+ (ITbl 0 0 EmptyOk (fromAssocs (BS1 0 (-1)) (BS1 (2^n-1) (Boundary $ n-1)) (999999) []))+ (ITbl 1 0 EmptyOk (fromAssocs Unit Unit (999999) []))+ EdgeWithSet+ Singleton+{-# NoInline forwardMinDist1 #-}++backtrackMinDist1 :: ScaleFunction -> Landscape -> Z:.TS1L Double:.U Double -> [Text]+backtrackMinDist1 scaleFunction landscape (Z:.ts1:.u) = unId $ axiom b+ where !(Z:.bt1:.b) = gMinDist (aMinDist scaleFunction landscape <|| aPretty scaleFunction landscape)+ (toBacktrack ts1 (undefined :: Id a -> Id a))+ (toBacktrack u (undefined :: Id a -> Id a))+ EdgeWithSet+ Singleton+ :: Z:.BT1L Double Text:.BTU Double Text+{-# NoInline backtrackMinDist1 #-}++-- | Count the number of co-optimals++minDistCount :: ScaleFunction -> Landscape -> Z:.TS1L (Double,Int):.U (Double,Int)+minDistCount scaleFunction landscape =+ let n = mutationCount landscape+ in mutateTablesST $ gMinDist (aMinDistCount scaleFunction landscape)+ (ITbl 0 0 EmptyOk (fromAssocs (BS1 0 (-1)) (BS1 (2^n-1) (Boundary $ n-1)) (999999,0) []))+ (ITbl 1 0 EmptyOk (fromAssocs Unit Unit (999999,0) []))+ EdgeWithSet+ Singleton+{-# NoInline minDistCount #-}++countBackMinDist1 :: ScaleFunction -> Landscape -> Z:.TS1L Double:.U Double -> [Integer]+countBackMinDist1 scaleFunction landscape (Z:.ts1:.u) = unId $ axiom b+ where !(Z:.bt1:.b) = gMinDist (aMinDist scaleFunction landscape <|| aCount landscape)+ (toBacktrack ts1 (undefined :: Id a -> Id a))+ (toBacktrack u (undefined :: Id a -> Id a))+ EdgeWithSet+ Singleton+ :: Z:.BT1L Double Integer:.BTU Double Integer+{-# NoInline countBackMinDist1 #-}++-- | Given the @Set1@ produced in @forwardMinDist1@ we can now extract the+-- co-optimal paths using the @Set1 -> ()@ index change.+--+-- TODO do we want this one explicitly or make life easy and just extract+-- from all @forwardMinDist1@ paths?++runCoOptDist :: ScaleFunction -> Landscape -> (Double,[Text])+runCoOptDist scaleFunction landscape = (unId $ axiom fwdu,bs)+ where !(Z:.fwd1:.fwdu) = forwardMinDist1 scaleFunction landscape+ bs = backtrackMinDist1 scaleFunction landscape (Z:.fwd1:.fwdu)+{-# NoInline runCoOptDist #-}++runCount :: ScaleFunction -> Landscape -> (Double,Int)+runCount scaleFunction landscape = (unId $ axiom fwdu)+ where !(Z:.fwd1:.fwdu) = minDistCount scaleFunction landscape+{-# NoInline runCount #-}++-- | Extract the individual partition scores.++boundaryPartFunFirst :: Maybe Int -> ScaleFunction -> Landscape -> [(Boundary First I,Log Double)]+boundaryPartFunFirst restrictStartNode scaleFunction landscape =+ let n = mutationCount landscape+ (Z:.sM:.bM) = mutateTablesST $ gMinDist (aInside restrictStartNode scaleFunction landscape)+ (ITbl 0 0 EmptyOk (fromAssocs (BS1 0 (-1)) (BS1 (2^n-1) (Boundary $ n-1)) (-999999) []))+ (ITbl 1 0 EmptyOk (fromAssocs (Boundary 0) (Boundary $ n-1) (-999999) []))+ EdgeWithSet+ Singleton+ :: Z:.TS1 (Log Double):.PF (Log Double)+ TW (ITbl _ _ _ pf) _ = bM+ bs' = assocs pf+ pssum = Numeric.Log.sum $ Prelude.map snd bs'+ bs = Prelude.map (second (/pssum)) bs'+ in bs+{-# NoInline boundaryPartFunFirst #-}++boundaryPartFunLast :: Maybe Int -> ScaleFunction -> Landscape -> BoundaryPart -- [(Boundary Last I,Log Double)]+boundaryPartFunLast restrictStartNode scaleFunction landscape =+ let n = mutationCount landscape+ (Z:.sM:.bM) = mutateTablesST $ gMinDist (aInside restrictStartNode scaleFunction landscape)+ (ITbl 0 0 EmptyOk (fromAssocs (BS1 0 (-1)) (BS1 (2^n-1) (Boundary $ n-1)) (-999999) []))+ (ITbl 1 0 EmptyOk (fromAssocs (Boundary 0) (Boundary $ n-1) (-999999) []))+ EdgeWithSet+ Singleton+ :: Z:.TS1L (Log Double):.PFL (Log Double)+ TW (ITbl _ _ _ pf) _ = bM+-- bs' = assocs pf+-- pssum = Numeric.Log.sum $ Prelude.map snd bs'+-- bs = Prelude.map (second (/pssum)) bs'+ in boundaryPart $ assocs pf -- bs+{-# NoInline boundaryPartFunLast #-}++data BoundaryPart = BoundaryPart+ { bpNormalized :: [(Boundary Last I, Log Double)]+ , bpUnnormalized :: [(Boundary Last I, Log Double)]+ , bpTotal :: Log Double+ }+ deriving (Show,Eq)++boundaryPart ps = BoundaryPart+ { bpNormalized = Prelude.map (second (/pssum)) ps+ , bpUnnormalized = ps+ , bpTotal = pssum+ }+ where pssum = Numeric.Log.sum $ Prelude.map snd ps+
+ BioInf/MutationOrder/RNA.hs view
@@ -0,0 +1,282 @@++-- | Here we collect the necessary data structures for the RNAs to be+-- compared. This data is serialized to disk once calculated, since it is+-- most likely the part that takes longest.+--+-- TODO if the vienna wrapper allows, we should parallelize the+-- calculations.+--+-- TODO nice interface counting up?++module BioInf.MutationOrder.RNA where++import Data.Aeson as DA+import Data.Bits+import Codec.Compression.GZip (compress,decompress)+import Control.Arrow (second)+import Control.DeepSeq+import Control.Parallel.Strategies+import Data.ByteString (ByteString)+import Data.Maybe (catMaybes)+import Data.Serialize+import Data.Vector.Serialize+import Data.Vector.Strategies+import Debug.Trace+import GHC.Generics+import qualified Data.ByteString.Char8 as BS+import qualified Data.ByteString.Lazy as BSL+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as VU+import System.IO.Unsafe (unsafePerformIO)+import qualified Data.HashMap.Strict as HM+import Data.Serialize.Instances+import Data.Text.Encoding (decodeUtf8, encodeUtf8)+import Data.Monoid+import Data.Char (isDigit)+import Data.Tuple (swap)++import qualified Data.Bijection.HashMap as B+import BioInf.ViennaRNA.Bindings+import qualified Data.PrimitiveArray as PA+import Biobase.Secondary.Diagrams (D1Secondary(..), mkD1S)++++-- | A single RNA with pre-calculated elements.+--+-- All calculations are done at 37 C.+--+-- TODO include the basepair probability matrix? Can we "compress" that+-- one?+--+-- We do not encode D1S into the json++data RNA = RNA+ { mutationSet :: !(VU.Vector (Int,Char))+ -- ^ we store just the mutation set, since this is more sparse and+ -- gives access to the mutational events.+ , primarySequence :: !ByteString+ -- ^ store RNA sequence too, for now+ , mfeStructure :: !ByteString+ -- ^ the mfe structure we get+ , mfeD1S :: !D1Secondary+ -- ^ efficient structure encoding+ , mfeEnergy :: !Double+ -- ^ mfe energy of the structure+ , centroidStructure :: !ByteString+ -- ^ the centroid structure+ , centroidD1S :: !D1Secondary+ -- ^ efficient centroid structure encoding+ , centroidEnergy :: !Double+ }+ deriving (Show,Eq,Generic)++instance NFData RNA+instance Serialize RNA++instance ToJSON RNA where+ toJSON RNA{..} =+ object [ "mutationSet" .= mutationSet+ , "primarySequence" .= decodeUtf8 primarySequence+ , "mfeStructure" .= decodeUtf8 mfeStructure+ , "mfeEnergy" .= mfeEnergy+ , "centroidStructure" .= decodeUtf8 centroidStructure+ , "centroidEnergy" .= centroidEnergy+ ]+ toEncoding RNA{..} =+ pairs ( "mutationSet" .= mutationSet+ <> "primarySequence" .= decodeUtf8 primarySequence+ <> "mfeStructure" .= decodeUtf8 mfeStructure+ <> "mfeEnergy" .= mfeEnergy+ <> "centroidStructure" .= decodeUtf8 centroidStructure+ <> "centroidEnergy" .= centroidEnergy+ )++instance FromJSON RNA where+ parseJSON (Object v) = do+ mutationSet <- v .: "mutationSet"+ primarySequence <- encodeUtf8 <$> v .: "primarySequence"+ let (e,s) = second BS.pack . unsafePerformIO . mfeTemp 37 $ BS.unpack primarySequence+ mfeStructure <- (fmap encodeUtf8 <$> v .:? "mfeStructure") .!= s+ mfeEnergy <- v .:? "mfeEnergy" .!= e+ let (ce,cs) = second BS.pack . unsafePerformIO . centroidTemp 37 $ BS.unpack primarySequence+ centroidStructure <- (fmap encodeUtf8 <$> v .:? "centroidStructure") .!= cs+ centroidEnergy <- v .:? "centroidEnergy" .!= ce+ let mfeD1S = bldD1S mfeStructure+ let centroidD1S = bldD1S centroidStructure+ return RNA{..}++bldD1S :: ByteString -> D1Secondary+bldD1S x = mkD1S (["()"::String], BS.unpack x)++-- | Given the primary sequence and the mutation set, fill the 'RNA'+-- structure.+--+-- NOTE This wraps some @ViennaRNA-bindings@ calls that are in @IO@.+--+-- TODO check if these calls are *really* thread-safe!++mkRNA+ :: Maybe (HM.HashMap ByteString QLine)+ -> ByteString+ -- ^ primary sequence of the *origin* RNA+ -> VU.Vector (Int,Char)+ -- ^ set of mutations compared to the origin+ -> RNA+mkRNA lkup inp' ms = RNA+ { mutationSet = ms+ , primarySequence = inp+ , mfeStructure = mS+ , mfeEnergy = mE+ , centroidStructure = cS+ , centroidEnergy = cE+ , mfeD1S = bldD1S mS+ , centroidD1S = bldD1S cS+ }+ where+ inp = insertMutations ms inp'+ ((mE,mS),(cE,cS)) = maybe calculateHere lookup lkup+ calculateHere = ( second BS.pack . unsafePerformIO . mfeTemp 37 $ BS.unpack inp+ , second BS.pack . unsafePerformIO . centroidTemp 37 $ BS.unpack inp+ )+ lookup lkup =+ case HM.lookup inp lkup of+ Nothing -> traceShow ("WARNING! have RNA lookup table but have to calculate!", inp) calculateHere+ Just QLine{..} -> (swap qlmfe,swap qlcentroid)++-- | Insert a set of mutations in a @ByteString@.++insertMutations :: VU.Vector (Int,Char) -> ByteString -> ByteString+insertMutations ms s' = VU.foldl' go s' ms+ where go s (k,c) =+ let (h,t) = BS.splitAt k s+ in BS.concat [h, BS.singleton c, BS.drop 1 t]++data Landscape = Landscape+ { rnas :: HM.HashMap (PA.BitSet PA.I) RNA+ -- ^ the individual RNA mutations. The index should be calculated from+ -- @linearIndex 0 high mutationSet@+ , mutationCount :: !Int+ -- ^ how many nucleotides are mutated in total+ , landscapeOrigin :: !ByteString+ -- ^ the ancestral sequence+ , landscapeDestination :: !ByteString+ -- ^ the final sequence+ , mutationPositions :: !(B.BimapHashMap Int Int)+ }+ deriving (Show,Eq,Generic)++instance NFData Landscape+instance Serialize Landscape++instance ToJSON Landscape where+ toJSON Landscape{..} =+ object [ "rnas" .= rnas+ , "mutationCount" .= mutationCount+ , "landscapeOrigin" .= decodeUtf8 landscapeOrigin+ , "landscapeDestination" .= decodeUtf8 landscapeDestination+ , "mutationPositions" .= mutationPositions+ ]+ toEncoding Landscape{..} =+ pairs ( "rnas" .= rnas+ <> "mutationCount" .= mutationCount+ <> "landscapeOrigin" .= decodeUtf8 landscapeOrigin+ <> "landscapeDestination" .= decodeUtf8 landscapeDestination+ <> "mutationPositions" .= mutationPositions+ )++instance FromJSON Landscape where+ parseJSON (Object v) = do+ rnas <- v .: "rnas"+ mutationCount <- v .: "mutationCount"+ landscapeOrigin <- encodeUtf8 <$> v .: "landscapeOrigin"+ landscapeDestination <- encodeUtf8 <$> v .: "landscapeDestination"+ mutationPositions <- v .: "mutationPositions"+ return Landscape{..}++-- |+--+-- TODO prime candidate for parallelization. ViennaRNA-bindings currently+-- does not allow parallel runs! It would be possible to consider+-- externalizing this, but for now we just run single-threaded.++createRNAlandscape :: Maybe (HM.HashMap ByteString QLine) -> Bool -> ByteString -> ByteString -> (Landscape, [(Int,ByteString)])+createRNAlandscape lkup verbose origin mutation = (ls, zipWith (\mm k -> (k,insertMutations mm origin)) mus [0..])+ where+ ls = Landscape+ { rnas = rs -- `using` (parVector chunkSize)+ , mutationCount = length . filter (>1) . map length $ pms+ , landscapeOrigin = origin+ , landscapeDestination = mutation+ , mutationPositions = mutbit+ }+ rs = HM.fromList . map pairWithBitSet $ zipWith talk mus [0..]+ talk s c = (if (c `mod` 1000 == 0 && verbose) then traceShow c else id) mkRNA lkup origin s+ mus = map (VU.fromList . catMaybes)+ . sequence+ $ pms+ -- possible mutations+ pms = zipWith3 genM (BS.unpack origin) (BS.unpack mutation) [0..]+ genM a b k | a==b = [Nothing]+ | otherwise = [Nothing,Just (k,b)]+ -- pair each @RNA@ with the correct bitset+ pairWithBitSet r = (calcBitSet zeroBits . map fst . VU.toList $ mutationSet r, r)+ -- calculate the bitset pattern for this mutation+ calcBitSet bs [] = bs+ calcBitSet bs (x':xs) =+ let x = maybe (error $ "calcBitSet") id $ B.lookupL mutbit x'+ in calcBitSet (bs `setBit` x) xs+ -- bijection between mutation position and bit position+ -- @BitSet Bit <-> Mutated Bit@+ mutbit = B.fromList+ . zipWith (flip (,)) [0 :: Int ..]+ . catMaybes $ zipWith3 genB (BS.unpack origin) (BS.unpack mutation) [0 :: Int ..]+ genB a b k | a == b = Nothing+ | otherwise = Just $ k++-- | Write a generated landscape to disk.++toFile :: FilePath -> Landscape -> IO ()+toFile fp = BSL.writeFile fp . compress . encodeLazy++toFileJSON :: FilePath -> Landscape -> IO ()+toFileJSON fp = BSL.writeFile fp . compress . DA.encode++fromFile :: FilePath -> IO Landscape+fromFile fp = (decodeLazy . decompress) <$> BSL.readFile fp >>= \case+ Left err -> error $ "BioInf.MutationOrder.RNA.fromFile: " ++ err+ Right ls -> return ls++fromFileJSON :: FilePath -> IO Landscape+fromFileJSON fp = (DA.eitherDecode' . decompress) <$> BSL.readFile fp >>= \case+ Left err -> error $ "BioInf.MutationOrder.RNA.fromFile: " ++ err+ Right ls -> return ls+++-- stupid parsing for quintuple rnafold lines++data QLine = QLine+ { qlSequence :: ByteString+ , qlmfe :: (ByteString,Double)+ , qlensemble :: (ByteString,Double)+ , qlcentroid :: (ByteString,Double)+ }+ deriving (Show)++++qlines f = do+ ls <- BS.lines <$> BS.readFile f+ return $ qlhm $ go ls+ where go [] = []+ go ls = let (hs,ts) = splitAt 5 ls+ in parseql hs : go ts+ parseql [s,m,e,c,_] =+ QLine s+ (stupid m)+ (stupid e)+ (stupid c)+ stupid bs = let (h:ts) = BS.words bs+ r = BS.dropWhile (\c -> not $ isDigit c || c=='-') $ BS.unwords ts+ in (h, read $ BS.unpack $ BS.takeWhile (\c -> isDigit c || c =='-' || c=='.') r)+ qlhm xs = HM.fromList $ map (\q -> (qlSequence q, q)) xs
+ LICENSE view
@@ -0,0 +1,675 @@+ GNU GENERAL PUBLIC LICENSE+ Version 3, 29 June 2007++ Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>+ Everyone is permitted to copy and distribute verbatim copies+ of this license document, but changing it is not allowed.++ Preamble++ The GNU General Public License is a free, copyleft license for+software and other kinds of works.++ The licenses for most software and other practical works are designed+to take away your freedom to share and change the works. By contrast,+the GNU General Public License is intended to guarantee your freedom to+share and change all versions of a program--to make sure it remains free+software for all its users. We, the Free Software Foundation, use the+GNU General Public License for most of our software; it applies also to+any other work released this way by its authors. You can apply it to+your programs, too.++ When we speak of free software, we are referring to freedom, not+price. Our General Public Licenses are designed to make sure that you+have the freedom to distribute copies of free software (and charge for+them if you wish), that you receive source code or can get it if you+want it, that you can change the software or use pieces of it in new+free programs, and that you know you can do these things.++ To protect your rights, we need to prevent others from denying you+these rights or asking you to surrender the rights. Therefore, you have+certain responsibilities if you distribute copies of the software, or if+you modify it: responsibilities to respect the freedom of others.++ For example, if you distribute copies of such a program, whether+gratis or for a fee, you must pass on to the recipients the same+freedoms that you received. You must make sure that they, too, receive+or can get the source code. And you must show them these terms so they+know their rights.++ Developers that use the GNU GPL protect your rights with two steps:+(1) assert copyright on the software, and (2) offer you this License+giving you legal permission to copy, distribute and/or modify it.++ For the developers' and authors' protection, the GPL clearly explains+that there is no warranty for this free software. For both users' and+authors' sake, the GPL requires that modified versions be marked as+changed, so that their problems will not be attributed erroneously to+authors of previous versions.++ Some devices are designed to deny users access to install or run+modified versions of the software inside them, although the manufacturer+can do so. This is fundamentally incompatible with the aim of+protecting users' freedom to change the software. The systematic+pattern of such abuse occurs in the area of products for individuals to+use, which is precisely where it is most unacceptable. Therefore, we+have designed this version of the GPL to prohibit the practice for those+products. If such problems arise substantially in other domains, we+stand ready to extend this provision to those domains in future versions+of the GPL, as needed to protect the freedom of users.++ Finally, every program is threatened constantly by software patents.+States should not allow patents to restrict development and use of+software on general-purpose computers, but in those that do, we wish to+avoid the special danger that patents applied to a free program could+make it effectively proprietary. To prevent this, the GPL assures that+patents cannot be used to render the program non-free.++ The precise terms and conditions for copying, distribution and+modification follow.++ TERMS AND CONDITIONS++ 0. Definitions.++ "This License" refers to version 3 of the GNU General Public License.++ "Copyright" also means copyright-like laws that apply to other kinds of+works, such as semiconductor masks.+ + "The Program" refers to any copyrightable work licensed under this+License. Each licensee is addressed as "you". "Licensees" and+"recipients" may be individuals or organizations.++ To "modify" a work means to copy from or adapt all or part of the work+in a fashion requiring copyright permission, other than the making of an+exact copy. The resulting work is called a "modified version" of the+earlier work or a work "based on" the earlier work.++ A "covered work" means either the unmodified Program or a work based+on the Program.++ To "propagate" a work means to do anything with it that, without+permission, would make you directly or secondarily liable for+infringement under applicable copyright law, except executing it on a+computer or modifying a private copy. Propagation includes copying,+distribution (with or without modification), making available to the+public, and in some countries other activities as well.++ To "convey" a work means any kind of propagation that enables other+parties to make or receive copies. Mere interaction with a user through+a computer network, with no transfer of a copy, is not conveying.++ An interactive user interface displays "Appropriate Legal Notices"+to the extent that it includes a convenient and prominently visible+feature that (1) displays an appropriate copyright notice, and (2)+tells the user that there is no warranty for the work (except to the+extent that warranties are provided), that licensees may convey the+work under this License, and how to view a copy of this License. If+the interface presents a list of user commands or options, such as a+menu, a prominent item in the list meets this criterion.++ 1. Source Code.++ The "source code" for a work means the preferred form of the work+for making modifications to it. "Object code" means any non-source+form of a work.++ A "Standard Interface" means an interface that either is an official+standard defined by a recognized standards body, or, in the case of+interfaces specified for a particular programming language, one that+is widely used among developers working in that language.++ The "System Libraries" of an executable work include anything, other+than the work as a whole, that (a) is included in the normal form of+packaging a Major Component, but which is not part of that Major+Component, and (b) serves only to enable use of the work with that+Major Component, or to implement a Standard Interface for which an+implementation is available to the public in source code form. A+"Major Component", in this context, means a major essential component+(kernel, window system, and so on) of the specific operating system+(if any) on which the executable work runs, or a compiler used to+produce the work, or an object code interpreter used to run it.++ The "Corresponding Source" for a work in object code form means all+the source code needed to generate, install, and (for an executable+work) run the object code and to modify the work, including scripts to+control those activities. However, it does not include the work's+System Libraries, or general-purpose tools or generally available free+programs which are used unmodified in performing those activities but+which are not part of the work. For example, Corresponding Source+includes interface definition files associated with source files for+the work, and the source code for shared libraries and dynamically+linked subprograms that the work is specifically designed to require,+such as by intimate data communication or control flow between those+subprograms and other parts of the work.++ The Corresponding Source need not include anything that users+can regenerate automatically from other parts of the Corresponding+Source.++ The Corresponding Source for a work in source code form is that+same work.++ 2. Basic Permissions.++ All rights granted under this License are granted for the term of+copyright on the Program, and are irrevocable provided the stated+conditions are met. This License explicitly affirms your unlimited+permission to run the unmodified Program. The output from running a+covered work is covered by this License only if the output, given its+content, constitutes a covered work. This License acknowledges your+rights of fair use or other equivalent, as provided by copyright law.++ You may make, run and propagate covered works that you do not+convey, without conditions so long as your license otherwise remains+in force. You may convey covered works to others for the sole purpose+of having them make modifications exclusively for you, or provide you+with facilities for running those works, provided that you comply with+the terms of this License in conveying all material for which you do+not control copyright. Those thus making or running the covered works+for you must do so exclusively on your behalf, under your direction+and control, on terms that prohibit them from making any copies of+your copyrighted material outside their relationship with you.++ Conveying under any other circumstances is permitted solely under+the conditions stated below. Sublicensing is not allowed; section 10+makes it unnecessary.++ 3. Protecting Users' Legal Rights From Anti-Circumvention Law.++ No covered work shall be deemed part of an effective technological+measure under any applicable law fulfilling obligations under article+11 of the WIPO copyright treaty adopted on 20 December 1996, or+similar laws prohibiting or restricting circumvention of such+measures.++ When you convey a covered work, you waive any legal power to forbid+circumvention of technological measures to the extent such circumvention+is effected by exercising rights under this License with respect to+the covered work, and you disclaim any intention to limit operation or+modification of the work as a means of enforcing, against the work's+users, your or third parties' legal rights to forbid circumvention of+technological measures.++ 4. Conveying Verbatim Copies.++ You may convey verbatim copies of the Program's source code as you+receive it, in any medium, provided that you conspicuously and+appropriately publish on each copy an appropriate copyright notice;+keep intact all notices stating that this License and any+non-permissive terms added in accord with section 7 apply to the code;+keep intact all notices of the absence of any warranty; and give all+recipients a copy of this License along with the Program.++ You may charge any price or no price for each copy that you convey,+and you may offer support or warranty protection for a fee.++ 5. Conveying Modified Source Versions.++ You may convey a work based on the Program, or the modifications to+produce it from the Program, in the form of source code under the+terms of section 4, provided that you also meet all of these conditions:++ a) The work must carry prominent notices stating that you modified+ it, and giving a relevant date.++ b) The work must carry prominent notices stating that it is+ released under this License and any conditions added under section+ 7. This requirement modifies the requirement in section 4 to+ "keep intact all notices".++ c) You must license the entire work, as a whole, under this+ License to anyone who comes into possession of a copy. This+ License will therefore apply, along with any applicable section 7+ additional terms, to the whole of the work, and all its parts,+ regardless of how they are packaged. This License gives no+ permission to license the work in any other way, but it does not+ invalidate such permission if you have separately received it.++ d) If the work has interactive user interfaces, each must display+ Appropriate Legal Notices; however, if the Program has interactive+ interfaces that do not display Appropriate Legal Notices, your+ work need not make them do so.++ A compilation of a covered work with other separate and independent+works, which are not by their nature extensions of the covered work,+and which are not combined with it such as to form a larger program,+in or on a volume of a storage or distribution medium, is called an+"aggregate" if the compilation and its resulting copyright are not+used to limit the access or legal rights of the compilation's users+beyond what the individual works permit. Inclusion of a covered work+in an aggregate does not cause this License to apply to the other+parts of the aggregate.++ 6. Conveying Non-Source Forms.++ You may convey a covered work in object code form under the terms+of sections 4 and 5, provided that you also convey the+machine-readable Corresponding Source under the terms of this License,+in one of these ways:++ a) Convey the object code in, or embodied in, a physical product+ (including a physical distribution medium), accompanied by the+ Corresponding Source fixed on a durable physical medium+ customarily used for software interchange.++ b) Convey the object code in, or embodied in, a physical product+ (including a physical distribution medium), accompanied by a+ written offer, valid for at least three years and valid for as+ long as you offer spare parts or customer support for that product+ model, to give anyone who possesses the object code either (1) a+ copy of the Corresponding Source for all the software in the+ product that is covered by this License, on a durable physical+ medium customarily used for software interchange, for a price no+ more than your reasonable cost of physically performing this+ conveying of source, or (2) access to copy the+ Corresponding Source from a network server at no charge.++ c) Convey individual copies of the object code with a copy of the+ written offer to provide the Corresponding Source. This+ alternative is allowed only occasionally and noncommercially, and+ only if you received the object code with such an offer, in accord+ with subsection 6b.++ d) Convey the object code by offering access from a designated+ place (gratis or for a charge), and offer equivalent access to the+ Corresponding Source in the same way through the same place at no+ further charge. You need not require recipients to copy the+ Corresponding Source along with the object code. If the place to+ copy the object code is a network server, the Corresponding Source+ may be on a different server (operated by you or a third party)+ that supports equivalent copying facilities, provided you maintain+ clear directions next to the object code saying where to find the+ Corresponding Source. Regardless of what server hosts the+ Corresponding Source, you remain obligated to ensure that it is+ available for as long as needed to satisfy these requirements.++ e) Convey the object code using peer-to-peer transmission, provided+ you inform other peers where the object code and Corresponding+ Source of the work are being offered to the general public at no+ charge under subsection 6d.++ A separable portion of the object code, whose source code is excluded+from the Corresponding Source as a System Library, need not be+included in conveying the object code work.++ A "User Product" is either (1) a "consumer product", which means any+tangible personal property which is normally used for personal, family,+or household purposes, or (2) anything designed or sold for incorporation+into a dwelling. In determining whether a product is a consumer product,+doubtful cases shall be resolved in favor of coverage. For a particular+product received by a particular user, "normally used" refers to a+typical or common use of that class of product, regardless of the status+of the particular user or of the way in which the particular user+actually uses, or expects or is expected to use, the product. A product+is a consumer product regardless of whether the product has substantial+commercial, industrial or non-consumer uses, unless such uses represent+the only significant mode of use of the product.++ "Installation Information" for a User Product means any methods,+procedures, authorization keys, or other information required to install+and execute modified versions of a covered work in that User Product from+a modified version of its Corresponding Source. The information must+suffice to ensure that the continued functioning of the modified object+code is in no case prevented or interfered with solely because+modification has been made.++ If you convey an object code work under this section in, or with, or+specifically for use in, a User Product, and the conveying occurs as+part of a transaction in which the right of possession and use of the+User Product is transferred to the recipient in perpetuity or for a+fixed term (regardless of how the transaction is characterized), the+Corresponding Source conveyed under this section must be accompanied+by the Installation Information. But this requirement does not apply+if neither you nor any third party retains the ability to install+modified object code on the User Product (for example, the work has+been installed in ROM).++ The requirement to provide Installation Information does not include a+requirement to continue to provide support service, warranty, or updates+for a work that has been modified or installed by the recipient, or for+the User Product in which it has been modified or installed. Access to a+network may be denied when the modification itself materially and+adversely affects the operation of the network or violates the rules and+protocols for communication across the network.++ Corresponding Source conveyed, and Installation Information provided,+in accord with this section must be in a format that is publicly+documented (and with an implementation available to the public in+source code form), and must require no special password or key for+unpacking, reading or copying.++ 7. Additional Terms.++ "Additional permissions" are terms that supplement the terms of this+License by making exceptions from one or more of its conditions.+Additional permissions that are applicable to the entire Program shall+be treated as though they were included in this License, to the extent+that they are valid under applicable law. If additional permissions+apply only to part of the Program, that part may be used separately+under those permissions, but the entire Program remains governed by+this License without regard to the additional permissions.++ When you convey a copy of a covered work, you may at your option+remove any additional permissions from that copy, or from any part of+it. (Additional permissions may be written to require their own+removal in certain cases when you modify the work.) You may place+additional permissions on material, added by you to a covered work,+for which you have or can give appropriate copyright permission.++ Notwithstanding any other provision of this License, for material you+add to a covered work, you may (if authorized by the copyright holders of+that material) supplement the terms of this License with terms:++ a) Disclaiming warranty or limiting liability differently from the+ terms of sections 15 and 16 of this License; or++ b) Requiring preservation of specified reasonable legal notices or+ author attributions in that material or in the Appropriate Legal+ Notices displayed by works containing it; or++ c) Prohibiting misrepresentation of the origin of that material, or+ requiring that modified versions of such material be marked in+ reasonable ways as different from the original version; or++ d) Limiting the use for publicity purposes of names of licensors or+ authors of the material; or++ e) Declining to grant rights under trademark law for use of some+ trade names, trademarks, or service marks; or++ f) Requiring indemnification of licensors and authors of that+ material by anyone who conveys the material (or modified versions of+ it) with contractual assumptions of liability to the recipient, for+ any liability that these contractual assumptions directly impose on+ those licensors and authors.++ All other non-permissive additional terms are considered "further+restrictions" within the meaning of section 10. If the Program as you+received it, or any part of it, contains a notice stating that it is+governed by this License along with a term that is a further+restriction, you may remove that term. If a license document contains+a further restriction but permits relicensing or conveying under this+License, you may add to a covered work material governed by the terms+of that license document, provided that the further restriction does+not survive such relicensing or conveying.++ If you add terms to a covered work in accord with this section, you+must place, in the relevant source files, a statement of the+additional terms that apply to those files, or a notice indicating+where to find the applicable terms.++ Additional terms, permissive or non-permissive, may be stated in the+form of a separately written license, or stated as exceptions;+the above requirements apply either way.++ 8. Termination.++ You may not propagate or modify a covered work except as expressly+provided under this License. Any attempt otherwise to propagate or+modify it is void, and will automatically terminate your rights under+this License (including any patent licenses granted under the third+paragraph of section 11).++ However, if you cease all violation of this License, then your+license from a particular copyright holder is reinstated (a)+provisionally, unless and until the copyright holder explicitly and+finally terminates your license, and (b) permanently, if the copyright+holder fails to notify you of the violation by some reasonable means+prior to 60 days after the cessation.++ Moreover, your license from a particular copyright holder is+reinstated permanently if the copyright holder notifies you of the+violation by some reasonable means, this is the first time you have+received notice of violation of this License (for any work) from that+copyright holder, and you cure the violation prior to 30 days after+your receipt of the notice.++ Termination of your rights under this section does not terminate the+licenses of parties who have received copies or rights from you under+this License. If your rights have been terminated and not permanently+reinstated, you do not qualify to receive new licenses for the same+material under section 10.++ 9. Acceptance Not Required for Having Copies.++ You are not required to accept this License in order to receive or+run a copy of the Program. Ancillary propagation of a covered work+occurring solely as a consequence of using peer-to-peer transmission+to receive a copy likewise does not require acceptance. However,+nothing other than this License grants you permission to propagate or+modify any covered work. These actions infringe copyright if you do+not accept this License. Therefore, by modifying or propagating a+covered work, you indicate your acceptance of this License to do so.++ 10. Automatic Licensing of Downstream Recipients.++ Each time you convey a covered work, the recipient automatically+receives a license from the original licensors, to run, modify and+propagate that work, subject to this License. You are not responsible+for enforcing compliance by third parties with this License.++ An "entity transaction" is a transaction transferring control of an+organization, or substantially all assets of one, or subdividing an+organization, or merging organizations. If propagation of a covered+work results from an entity transaction, each party to that+transaction who receives a copy of the work also receives whatever+licenses to the work the party's predecessor in interest had or could+give under the previous paragraph, plus a right to possession of the+Corresponding Source of the work from the predecessor in interest, if+the predecessor has it or can get it with reasonable efforts.++ You may not impose any further restrictions on the exercise of the+rights granted or affirmed under this License. For example, you may+not impose a license fee, royalty, or other charge for exercise of+rights granted under this License, and you may not initiate litigation+(including a cross-claim or counterclaim in a lawsuit) alleging that+any patent claim is infringed by making, using, selling, offering for+sale, or importing the Program or any portion of it.++ 11. Patents.++ A "contributor" is a copyright holder who authorizes use under this+License of the Program or a work on which the Program is based. The+work thus licensed is called the contributor's "contributor version".++ A contributor's "essential patent claims" are all patent claims+owned or controlled by the contributor, whether already acquired or+hereafter acquired, that would be infringed by some manner, permitted+by this License, of making, using, or selling its contributor version,+but do not include claims that would be infringed only as a+consequence of further modification of the contributor version. For+purposes of this definition, "control" includes the right to grant+patent sublicenses in a manner consistent with the requirements of+this License.++ Each contributor grants you a non-exclusive, worldwide, royalty-free+patent license under the contributor's essential patent claims, to+make, use, sell, offer for sale, import and otherwise run, modify and+propagate the contents of its contributor version.++ In the following three paragraphs, a "patent license" is any express+agreement or commitment, however denominated, not to enforce a patent+(such as an express permission to practice a patent or covenant not to+sue for patent infringement). To "grant" such a patent license to a+party means to make such an agreement or commitment not to enforce a+patent against the party.++ If you convey a covered work, knowingly relying on a patent license,+and the Corresponding Source of the work is not available for anyone+to copy, free of charge and under the terms of this License, through a+publicly available network server or other readily accessible means,+then you must either (1) cause the Corresponding Source to be so+available, or (2) arrange to deprive yourself of the benefit of the+patent license for this particular work, or (3) arrange, in a manner+consistent with the requirements of this License, to extend the patent+license to downstream recipients. "Knowingly relying" means you have+actual knowledge that, but for the patent license, your conveying the+covered work in a country, or your recipient's use of the covered work+in a country, would infringe one or more identifiable patents in that+country that you have reason to believe are valid.+ + If, pursuant to or in connection with a single transaction or+arrangement, you convey, or propagate by procuring conveyance of, a+covered work, and grant a patent license to some of the parties+receiving the covered work authorizing them to use, propagate, modify+or convey a specific copy of the covered work, then the patent license+you grant is automatically extended to all recipients of the covered+work and works based on it.++ A patent license is "discriminatory" if it does not include within+the scope of its coverage, prohibits the exercise of, or is+conditioned on the non-exercise of one or more of the rights that are+specifically granted under this License. You may not convey a covered+work if you are a party to an arrangement with a third party that is+in the business of distributing software, under which you make payment+to the third party based on the extent of your activity of conveying+the work, and under which the third party grants, to any of the+parties who would receive the covered work from you, a discriminatory+patent license (a) in connection with copies of the covered work+conveyed by you (or copies made from those copies), or (b) primarily+for and in connection with specific products or compilations that+contain the covered work, unless you entered into that arrangement,+or that patent license was granted, prior to 28 March 2007.++ Nothing in this License shall be construed as excluding or limiting+any implied license or other defenses to infringement that may+otherwise be available to you under applicable patent law.++ 12. No Surrender of Others' Freedom.++ If conditions are imposed on you (whether by court order, agreement or+otherwise) that contradict the conditions of this License, they do not+excuse you from the conditions of this License. If you cannot convey a+covered work so as to satisfy simultaneously your obligations under this+License and any other pertinent obligations, then as a consequence you may+not convey it at all. For example, if you agree to terms that obligate you+to collect a royalty for further conveying from those to whom you convey+the Program, the only way you could satisfy both those terms and this+License would be to refrain entirely from conveying the Program.++ 13. Use with the GNU Affero General Public License.++ Notwithstanding any other provision of this License, you have+permission to link or combine any covered work with a work licensed+under version 3 of the GNU Affero General Public License into a single+combined work, and to convey the resulting work. The terms of this+License will continue to apply to the part which is the covered work,+but the special requirements of the GNU Affero General Public License,+section 13, concerning interaction through a network will apply to the+combination as such.++ 14. Revised Versions of this License.++ The Free Software Foundation may publish revised and/or new versions of+the GNU General Public License from time to time. Such new versions will+be similar in spirit to the present version, but may differ in detail to+address new problems or concerns.++ Each version is given a distinguishing version number. If the+Program specifies that a certain numbered version of the GNU General+Public License "or any later version" applies to it, you have the+option of following the terms and conditions either of that numbered+version or of any later version published by the Free Software+Foundation. If the Program does not specify a version number of the+GNU General Public License, you may choose any version ever published+by the Free Software Foundation.++ If the Program specifies that a proxy can decide which future+versions of the GNU General Public License can be used, that proxy's+public statement of acceptance of a version permanently authorizes you+to choose that version for the Program.++ Later license versions may give you additional or different+permissions. However, no additional obligations are imposed on any+author or copyright holder as a result of your choosing to follow a+later version.++ 15. Disclaimer of Warranty.++ THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY+APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT+HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY+OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO,+THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR+PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM+IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF+ALL NECESSARY SERVICING, REPAIR OR CORRECTION.++ 16. Limitation of Liability.++ IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS+THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY+GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE+USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF+DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD+PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS),+EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF+SUCH DAMAGES.++ 17. Interpretation of Sections 15 and 16.++ If the disclaimer of warranty and limitation of liability provided+above cannot be given local legal effect according to their terms,+reviewing courts shall apply local law that most closely approximates+an absolute waiver of all civil liability in connection with the+Program, unless a warranty or assumption of liability accompanies a+copy of the Program in return for a fee.++ END OF TERMS AND CONDITIONS++ How to Apply These Terms to Your New Programs++ If you develop a new program, and you want it to be of the greatest+possible use to the public, the best way to achieve this is to make it+free software which everyone can redistribute and change under these terms.++ To do so, attach the following notices to the program. It is safest+to attach them to the start of each source file to most effectively+state the exclusion of warranty; and each file should have at least+the "copyright" line and a pointer to where the full notice is found.++ <one line to give the program's name and a brief idea of what it does.>+ Copyright (C) <year> <name of author>++ This program is free software: you can redistribute it and/or modify+ it under the terms of the GNU General Public License as published by+ the Free Software Foundation, either version 3 of the License, or+ (at your option) any later version.++ This program is distributed in the hope that it will be useful,+ but WITHOUT ANY WARRANTY; without even the implied warranty of+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the+ GNU General Public License for more details.++ You should have received a copy of the GNU General Public License+ along with this program. If not, see <http://www.gnu.org/licenses/>.++Also add information on how to contact you by electronic and paper mail.++ If the program does terminal interaction, make it output a short+notice like this when it starts in an interactive mode:++ <program> Copyright (C) <year> <name of author>+ This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'.+ This is free software, and you are welcome to redistribute it+ under certain conditions; type `show c' for details.++The hypothetical commands `show w' and `show c' should show the appropriate+parts of the General Public License. Of course, your program's commands+might be different; for a GUI interface, you would use an "about box".++ You should also get your employer (if you work as a programmer) or school,+if any, to sign a "copyright disclaimer" for the program, if necessary.+For more information on this, and how to apply and follow the GNU GPL, see+<http://www.gnu.org/licenses/>.++ The GNU General Public License does not permit incorporating your program+into proprietary programs. If your program is a subroutine library, you+may consider it more useful to permit linking proprietary applications with+the library. If this is what you want to do, use the GNU Lesser General+Public License instead of this License. But first, please read+<http://www.gnu.org/philosophy/why-not-lgpl.html>.+
+ MutationOrder.cabal view
@@ -0,0 +1,166 @@+name: MutationOrder+version: 0.0.0.1+author: Maria Beatriz Walter Costa, Christian Hoener zu Siederdissen, 2017+copyright: Christian Hoener zu Siederdissen, 2017+homepage: https://github.com/choener/MutationOrder+bug-reports: https://github.com/choener/MutationOrder/issues+maintainer: choener@bioinf.uni-leipzig.de+category: Bioinformatics+license: GPL-3+license-file: LICENSE+build-type: Simple+stability: experimental+cabal-version: >= 1.10.0+tested-with: GHC == 7.10.3, GHC == 8.0.1+synopsis: Most likely order of mutation events in RNA+description:+ Determine the most likely order in which single nucleotide+ mutations happened between two RNA sequences.+ .+ Developed to analyse the @HAR 1@ region.+ .+ As long as the two input RNAs are small enough enough (couple+ hundred nucleotides) and the number of mutations is small+ enough (around 20-26, since the algorithm is exponential in+ this number) the algorithm should work for similar problems+ without changes.++++Extra-Source-Files:+ README.md+ changelog.md++++flag debug+ description: Enable bounds checking and various other debug operations at the cost of a significant performance penalty.+ default: False+ manual: True++flag debugoutput+ description: Enable debug output, which spams the screen full of index information+ default: False+ manual: True++++library+ build-depends: base >= 4.7 && < 5.0+ , aeson >= 1.1+ , bytestring+ , cereal >= 0.5+ , cereal-vector >= 0.2+ , containers+ , deepseq >= 1.4+ , directory+ , filepath+ , log-domain >= 0.10+ , parallel >= 3.2+ , serialize-instances >= 0.1+ , text >= 1.0+ , unordered-containers >= 0.2.7+ , vector >= 0.11+ , vector-strategies >= 0.4+ , zlib >= 0.6+ --+ , ADPfusion == 0.5.2.*+ , ADPfusionSet == 0.0.0.*+ , bimaps == 0.1.0.*+ , BiobaseXNA == 0.9.3.*+ , DPutils == 0.0.1.*+ , FormalGrammars == 0.3.1.*+ , PrimitiveArray == 0.8.0.*+ , PrimitiveArray-Pretty == 0.0.0.*+ , ShortestPathProblems == 0.0.0.*+ , ViennaRNA-bindings == 0.233.1.*+ exposed-modules:+ BioInf.MutationOrder+ BioInf.MutationOrder.EdgeProb+ BioInf.MutationOrder.MinDist+ BioInf.MutationOrder.RNA+ default-extensions: BangPatterns+ , CPP+ , DeriveDataTypeable+ , DeriveGeneric+ , FlexibleContexts+ , GADTs+ , LambdaCase+ , MultiParamTypeClasses+ , OverloadedStrings+ , QuasiQuotes+ , RecordWildCards+ , ScopedTypeVariables+ , TemplateHaskell+ , TupleSections+ , TypeFamilies+ , TypeOperators+ default-language:+ Haskell2010+ ghc-options:+ -O2 -funbox-strict-fields+ if flag(debug)+ cpp-options: -DADPFUSION_CHECKS+ ghc-options: -fno-ignore-asserts -O0+ if flag(debugoutput)+ cpp-options: -DADPFUSION_DEBUGOUTPUT++++executable MutationOrder+ build-depends: base+ , bytestring+ , cmdargs >= 0.10+ , filepath+ --+ , MutationOrder+ hs-source-dirs:+ src+ default-extensions: BangPatterns+ , DeriveDataTypeable+ , RecordWildCards+ main-is:+ MutationOrder.hs+ default-language:+ Haskell2010+ ghc-options:+ -O2 -rtsopts+-- -threaded++++test-suite properties+ type:+ exitcode-stdio-1.0+ main-is:+ properties.hs+ ghc-options:+ -threaded -rtsopts -with-rtsopts=-N+ hs-source-dirs:+ tests+ default-language:+ Haskell2010+ default-extensions: BangPatterns+ , CPP+ , FlexibleContexts+ , FlexibleInstances+ , MultiParamTypeClasses+ , ScopedTypeVariables+ , TemplateHaskell+ , TypeFamilies+ , TypeOperators+ , TypeSynonymInstances+ build-depends: base+ , QuickCheck+ , tasty >= 0.11+ , tasty-quickcheck >= 0.8+ , tasty-th >= 0.1+ , vector+ --+ , MutationOrder+++source-repository head+ type: git+ location: git://github.com/choener/MutationOrder+
+ README.md view
@@ -0,0 +1,48 @@+[](https://travis-ci.org/choener/MutationOrder)++Determine the most likely order of mutations from one RNA sequence to another.++1. Walter Costa, Maria Beatriz and Hoener zu Siederdissen, Christian and Tulpan, Dan and Stadler, Peter F. and Nowick, Katja + *Uncovering the Structural Evolution of the Human Accelerated Region 1* + 2017, submitted + [preprint](http://www.bioinf.uni-leipzig.de/~choener/pdfs/wal-hoe-2017.pdf) ++# Usage instructions++We assume that you have two Fasta files, *from.fa* and *to.fa* but they can be+named however is convenient (say *chimp.fa* and *human.fa*). Each file has to+contain exactly one sequence and both sequences have to be of the same length.++We then run++```./MutationOrder --workdb from-to.json.gz --scoretype pairdistcen --onlypositive --outputprefix test```++This will generate ```test.run```, ```test-edge.eps```, and+```test-meaorder.eps```. The ```test.run``` file provides extensive output of+the optimal path, the first-last probabilities, the edge probabilities, and the+mea output. The two ```eps``` files give a graphical representation of the edge+probabilities, for the ```meaorder``` in order of the path of maximum expected+accuracy.++The ```--scoretype``` allows for ```mfe```, ```centroid```, ```pairdistcen```,+and ```pairdistmfe```, which analyse possible evoluationary paths according to+mfe energy, centroid energy, smallest base pair distances for each step in the+```cen```troid or ```mfe``` case.++++# Installation++Follow [this+link](http://www.bioinf.uni-leipzig.de/~choener/software/MutationOrder.html) to+the bottom of the page. Binaries are available for download and installation+from sources via *Haskell Stack* are described.+++#### Contact++Christian Hoener zu Siederdissen +Leipzig University, Leipzig, Germany +choener@bioinf.uni-leipzig.de +http://www.bioinf.uni-leipzig.de/~choener/ +
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ changelog.md view
@@ -0,0 +1,6 @@+0.0.0.1+-------++- initial checkin+- travis-ci integration+
+ src/MutationOrder.hs view
@@ -0,0 +1,100 @@++{-# Options_GHC -fno-cse #-}++module Main where++import System.Console.CmdArgs+import System.FilePath+import qualified Data.ByteString.Char8 as BS+import Control.Monad++import BioInf.MutationOrder+import BioInf.MutationOrder.RNA (createRNAlandscape)++data ScoreType+ = Mfe+ | Centroid+ | PairDistMfe+ | PairDistCen+ deriving (Show,Data,Typeable)++data Options+ = Options+ { infiles :: [FilePath]+ , workdb :: FilePath+ , temperature :: Double+ , fillweight :: FillWeight+ , fillstyle :: FillStyle+ , cooptcount :: Int+ , cooptprint :: Int+ , outprefix :: FilePath+ , scoretype :: ScoreType+ , positivesquared :: Bool+ , onlypositive :: Bool+ , equalStart :: Bool+ , posscaled :: Maybe (Double,Double)+ , lkupfile :: Maybe FilePath+ }+ | GenSequences+ { infiles :: [FilePath]+ }+ deriving (Show,Data,Typeable)++oOptions = Options+ { infiles = def &= args+ , workdb = "work.db" &= help "name of the database to store intermediates in"+ , temperature = 1.0 &= help "lower temperatures favor the more optimal paths, defaults to 1.0"+ , fillweight = FWlog+ , fillstyle = FSfull+ , cooptcount = 100000+ , cooptprint = 2+ , outprefix = "tmp"+ , scoretype = Centroid &= help "choose 'mfe', 'centroid', 'pairdistmfe', or 'pairdistcen' for the evaluation of each mutational step"+ , positivesquared = False &= help "square positive energies to penalize worse structures"+ , onlypositive = False &= help "minimize only over penalties, not energy gains"+ , equalStart = False+ , posscaled = Nothing+ , lkupfile = Nothing+ }++oGenSequences = GenSequences+ { infiles = def &= args+ }++main :: IO ()+main = do+ o <- cmdArgs $ modes [oOptions, oGenSequences] &= verbosity+ case o of+ Options{} -> mainProgram o+ GenSequences{} -> genSequences o++genSequences o = do+ let GenSequences{..} = o+ ancestral <- stupidReader $ infiles !! 0+ current <- stupidReader $ infiles !! 1+ let ls = snd $ createRNAlandscape Nothing False ancestral current+ forM_ ls $ \(k,sq) -> BS.putStrLn sq+ return ()++mainProgram oOptions = do+ let Options{..} = oOptions+ isL <- isLoud+ let fwdScaleFunction+ = (if positivesquared then squaredPositive else id)+ . (maybe id (uncurry posScaled) posscaled)+ . (if onlypositive then (scaleByFunction (max 0)) else id)+ $ (case scoretype of Mfe -> mfeDelta+ Centroid -> centroidDelta+ PairDistMfe -> basepairDistanceMFE+ PairDistCen -> basepairDistanceCentroid)+ let insideScaleFunction+ = scaleTemperature temperature+ . (if positivesquared then squaredPositive else id)+ . (maybe id (uncurry posScaled) posscaled)+ . (if onlypositive then (scaleByFunction (max 0)) else id)+ $ (case scoretype of Mfe -> mfeDelta+ Centroid -> centroidDelta+ PairDistMfe -> basepairDistanceMFE+ PairDistCen -> basepairDistanceCentroid)+ runMutationOrder isL fillweight fillstyle fwdScaleFunction insideScaleFunction cooptcount cooptprint lkupfile outprefix workdb temperature equalStart infiles+
+ tests/properties.hs view
@@ -0,0 +1,14 @@++module Main where++import Test.Tasty+import Test.Tasty.TH++++main :: IO ()+main = do+ defaultMain $ testGroup ""+ [+ ]+