MutationOrder 0.0.0.2 → 0.0.1.0
raw patch · 9 files changed
+1168/−86 lines, 9 filesdep +OrderedBitsdep +attoparsecdep +bytestring-trie
Dependencies added: OrderedBits, attoparsec, bytestring-trie, errors, filemanip, hashable, lens, mtl, split
Files
- BioInf/MutationOrder.hs +319/−16
- BioInf/MutationOrder/BackMutations.hs +390/−0
- BioInf/MutationOrder/MinDist.hs +25/−25
- BioInf/MutationOrder/RNA.hs +84/−8
- BioInf/MutationOrder/SequenceDB.hs +167/−0
- MutationOrder.cabal +30/−4
- README.md +22/−2
- changelog.md +6/−0
- src/MutationOrder.hs +125/−31
BioInf/MutationOrder.hs view
@@ -29,57 +29,109 @@ , ScaleFunction (..) ) where -import qualified Data.Vector.Unboxed as VU-import Data.Tuple (swap)-import Control.Monad (unless,forM_,when)+import Control.Arrow (first,second)+import Control.Error+import Control.Monad.IO.Class (liftIO)+import Control.Monad (unless,forM_,when,forM)+import Control.Lens+import Data.ByteString.Strict.Lens import Data.Bits import Data.ByteString (ByteString)+import Data.Char (toUpper)+import Data.Data import Data.Function (on)-import Data.List (groupBy,sortBy)+import Data.List (groupBy,sortBy,foldl',(\\),sort)+import Data.List.Split (chunksOf) import Data.Ord (comparing)+import Data.Tuple (swap)+import Debug.Trace import Numeric.Log+import qualified Control.Parallel.Strategies as Par 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 qualified Data.Trie as Trie+import qualified Data.Vector.Unboxed as VU 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 System.IO (withFile,IOMode(WriteMode),hPutStrLn,Handle,stderr)+import Text.Printf import ADP.Fusion.Term.Edge.Type (From(..),To(..))+import Biobase.Secondary.Diagrams (d1Distance) 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 Data.PrimitiveArray as PA import qualified ShortestPath.SHP.Edge.MinDist as SHP-import Biobase.Secondary.Diagrams (d1Distance)+import Data.Bits.Ordered import BioInf.MutationOrder.EdgeProb import BioInf.MutationOrder.MinDist import BioInf.MutationOrder.RNA+import BioInf.MutationOrder.SequenceDB+import qualified BioInf.MutationOrder.BackMutations as BM -runMutationOrder verbose fw fs fwdScaleFunction probScaleFunction cooptCount cooptPrint lkupFile outprefix workdb temperature equalStart [ancestralFP,currentFP] = do+data ScoreType+ = Mfe+ | Centroid+ | PairDistMfe+ | PairDistCen+ | PairDistMfeTgt+ | PairDistCenTgt+ deriving (Show,Data,Typeable)++++runMutationOrder verbose fw fs scoretype positivesquared posscaled onlypositive cooptCount cooptPrint lkupFile outprefix workdb temperature equalStart [ancestralFP,currentFP] everyKth = 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+ hPrintf stderr "%s.run exists, ending now!\n" outprefix exitSuccess withFile (outprefix ++ ".run") WriteMode $ \oH -> do- printf "%s.run job started!\n" outprefix+ hPrintf stderr "%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+ hPrintf stderr "prepare to load dump file\n"+ !ls <- withDumpFile oH workdb ancestral current . fst $ createRNAlandscape lkup verbose ancestral current+ hPrintf stderr "dump file loaded\n"+ -- final state in the rna landscape+ let !tgtRNA = rnas ls HM.! (2^ mutationCount ls - 1)+ hPrintf stderr "rna db loaded, tgt rna set\n"+ 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+ PairDistMfeTgt → bpMFEDistToExtant tgtRNA+ PairDistCenTgt → bpCentroidDistToExtant tgtRNA+ )+ let probScaleFunction+ = 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+ PairDistMfeTgt → bpMFEDistToExtant tgtRNA+ PairDistCenTgt → bpCentroidDistToExtant tgtRNA+ ) let mpks = sortBy (comparing snd) . B.toList $ mutationPositions ls let bitToNuc = M.fromList $ map (swap . first (+1)) mpks let nn = length mpks@@ -88,14 +140,23 @@ -- -- Run co-optimal lowest energy changes --- let (e,bs) = runCoOptDist fwdScaleFunction ls+ hPrintf stderr "starting forward phase\n"+ let (!e,bs) = runCoOptDist fwdScaleFunction ls+ hPrintf stderr "forward phase done with %f\n" e+-- forM_ bs $ \b → hPrintf stderr "%s\n" (show b) 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 stderr "starting coopthisto\n"+ let takeKth k [] = []+ takeKth k (x:xs) = x : takeKth k xs+ let ch = coopthisto $ map snd $ takeKth everyKth $ take cooptCount $ printbs ++ countbs+ hPrintf stderr (show ch ++ "\n")+ -- TODO here we can now do a histogram with printbs and countbs 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)+ forM_ printbs (T.hPutStrLn oH . fst) hPrintf oH "%s\n\n" $ replicate 80 '=' -- -- Run @First@ probability algorithm to determine the probability for@@ -156,6 +217,7 @@ -- colMarginals gives the total probability that the mutation order -- ends with this mutation. let colMarginals = M.mapKeysWith (+) snd firstlastP+ hPrintf oH "lnZ: %8.5f\n" $ ln firstlastZ hPrintf oH " " forM_ (M.elems bitToNuc) $ \mut -> hPrintf oH "%6d " mut hPrintf oH " Σ\n"@@ -290,6 +352,24 @@ (SHP.BTedge (From ff:.To tt),mut) -> prettyPrint mut (Just tt) hPutStrLn oH "" hPutStrLn oH ""+ -- NEW: coopt histograms+ hPutStrLn oH "Co-optimality histograms"+ hPrintf oH " "+ forM_ [1.. mutationCount ls] (hPrintf oH " %10d")+ hPrintf oH "\n"+ forM_ ch $ \(h,as,rs) → do+ hPrintf oH "%5d" $ h+1+ -- absolutes+ VU.mapM_ (\c → hPrintf oH " %10d" c) as+ -- relatives+ hPrintf oH "\n r"+ VU.mapM_ (\c → hPrintf oH " %10.8f" c) rs+-- -- log-scale (1 + relative count)+-- hPrintf oH "\n lnR "+-- VU.mapM_ (\c → hPrintf oH " %10.8f" (log $ 1+c)) rs+ hPrintf oH "\n"+ hPrintf oH "\n"+-- -- MEA order let meaOrder = let go = \case SHP.BTnode (_:.To n) -> n SHP.BTedge (From ff:.To tt) -> tt@@ -313,6 +393,31 @@ -} {-# NoInline runMutationOrder #-} +-- | Histogram that provides for each mutation position the relative frequency+-- with which the mutation shows up as the k'th one in the order.++coopthisto+ ∷ [[Int]]+ -- ^ list of permutation orders @[ [116, 10, 48, 30] , [10, 30, 48, 116], ... ]@+ → [(Int,VU.Vector Int, VU.Vector Double)]+ -- ^ list of mutation, absolute frequency of position, relative frequency of position+coopthisto as@(a':_) = es -- traceShow (keys,as,bs,cs,ds) $ es+ where+ -- just so that we know+ keys = sort a'+ -- add positional information: each permutation order has @(,index)@ added.+ -- then we can concatenate all lists and get @[ (116,1), (10,2), ... (116,4), ...]@+ rks = reverse [1 .. length keys]+ bs = concatMap (\a → zip a rks) as+ -- into hashmap, (mutation, order position) sum up how often seen+ cs = HM.fromListWith (+) $ map (,1∷Int) bs+ -- absolute frequencies+ ds = [ (k, VU.fromList [ HM.lookupDefault 0 (k, i) cs+ | i ← [1 .. length keys] ])+ | k ← keys ]+ -- attach relative frequencies+ es = [ (k, abs, VU.map (\q → fromIntegral q / s) abs) | (k,abs) ← ds, let (s∷Double) = fromIntegral $ VU.sum abs ]+ posScaled :: Double -> Double -> ScaleFunction -> ScaleFunction posScaled l s = scaleByFunction go where go d | d >= l = d ** s@@ -320,7 +425,7 @@ {-# Inline go #-} {-# Inlinable posScaled #-} --- | Basepair distance+-- ** Basepair distance, between neighbors basepairDistanceMFE :: ScaleFunction basepairDistanceMFE frna trna = fromIntegral $ d1Distance (mfeD1S frna) (mfeD1S trna)@@ -328,6 +433,16 @@ basepairDistanceCentroid :: ScaleFunction basepairDistanceCentroid frna trna = fromIntegral $ d1Distance (centroidD1S frna) (centroidD1S trna) +-- ** Basepair distance, from currently newest to extant++bpMFEDistToExtant ∷ RNA → ScaleFunction+bpMFEDistToExtant extant frna trna = fromIntegral $ d1Distance e (mfeD1S trna) - d1Distance e (mfeD1S frna)+ where e = mfeD1S extant++bpCentroidDistToExtant ∷ RNA → ScaleFunction+bpCentroidDistToExtant extant frna trna = fromIntegral $ d1Distance e (mfeD1S trna) - d1Distance e (mfeD1S frna)+ where e = mfeD1S extant+ -- | Scale function for normal mfe delta energies mfeDelta :: ScaleFunction@@ -408,4 +523,192 @@ hPrintf oH "database %s does not exist! Folding all intermediate structures. This may take a while!\n" fp toFileJSON fp l return l++-- | This function will run a subset of the possible backmutations. In+-- particular, only those mutations for one particular backmutation column are+-- being used.+--+-- TODO less explicit transformer stack!++runBackmutationVariants+ ∷ BM.ScaleFunction Double+ → FilePath+ -- ^ where the work db lives+ → [Char]+ -- ^+ → Ancestral+ -- ^ The ancestral sequence from which we mutate away+ → Extant+ -- ^ The sequence to which to mutate to+ → Int+ -- ^ The backmutation / intermediate mutation we look at. Indexed with @[1..sequence length]@.+ -- TODO wrap in newtype that enforces this. We have some index structure saying "Start at 1".+ → ExceptT String IO ()+runBackmutationVariants scaleFun workdb alphabet ancestral extant ipos' = do+ let ipos = ipos' - 1+ -- error guarding+ unless (ipos >=0) $ throwE "ipos can not be negative"+ unless (ipos < (BS.length $ getAncestral ancestral)) $ throwE "ipos larger than sequence length"+ unless (BS.length (getAncestral ancestral) == BS.length (getExtant extant)) $ throwE "ancestral and extant sequence do not have equal length"+ -- Load all sequences for the original problem -- they are needed anyway+ -- TODO i think @origSeqs'@ does not hold the original sequences, need to check tonight+ (seqCount, origSeqs', variants) ← createRNAlandscape2 alphabet (Right [ipos]) [] ancestral extant+ let origSeqs = Trie.fromList [ (s,()) | s ← origSeqs' ]+ varSeqs ← forM (alphabet \\ [getAncestral ancestral `BS.index` ipos, getExtant extant `BS.index` ipos]) $ \v → do+ let vsqs = Trie.fromList [ (s,()) | (i,c,ss) ← variants, c == v, s ← ss ]+ return (v, vsqs)+ -- liftIO $ print workdb+ -- TODO do only a single pass over the data+ -- origStrs ← liftIO $ filter (\r → rnaFoldSequence r `Trie.member` origSeqs) <$> map rna2dna <$> readRNAfoldFiles workdb+ allStrs ← liftIO $ filter (\r → rnaFoldSequence r `Trie.member` origSeqs+ || or [rnaFoldSequence r `Trie.member` vs | (_,vs) ← varSeqs] )+ <$> map rna2dna <$> readRNAfoldFiles workdb+ let (rnas,_,_,_) = genSet ancestral extant Nothing allStrs+ -- The @ipos@ declares how many variants we have+ forM_ (alphabet \\ [getAncestral ancestral `BS.index` ipos, getExtant extant `BS.index` ipos]) $ \v → do+ let varSeqs = Trie.fromList [ (s,()) | (i,c,ss) ← variants, c == v, s ← ss ]+ varStrs ← liftIO $ filter (\r → rnaFoldSequence r `Trie.member` varSeqs) <$> map rna2dna <$> readRNAfoldFiles workdb+ let (ntrs,iposbitset,numMuts,mutpos) = genSet ancestral extant (Just (ipos,v)) allStrs+ --liftIO $ print (numMuts, iposbitset, rnas, ntrs, HM.size rnas, HM.size ntrs)+ let fmd = BM.forwardMinDist numMuts+ scaleFun+ iposbitset+ rnas+ ntrs+ let bts = BM.backtrackMinDist1 numMuts+ scaleFun+ iposbitset+ ipos'+ rnas+ ntrs+ mutpos+ fmd+ let evi = BM.forwardEvidence numMuts+ (partfun' scaleFun)+ iposbitset+ rnas+ ntrs+ liftIO $ printf "Unobserved Mutation Type: %s Position: %3d Nucleotide: %c Delta: %5.1f lnZ: %10.2f\n"+ (if iposbitset < 0 then "BM" else "IM" :: String)+ ipos' v+ (BM.forwardMinDistValue fmd)+ (ln $ BM.forwardEvidenceValue evi)+ liftIO $ mapM_ T.putStrLn $ take 1 bts+ return ()+ return ()++mfeDelta' :: Bool → Bool → BM.ScaleFunction Double+mfeDelta' mx sq frna trna = guardPositives mx sq $ t - f+ where t = rnaFoldMFEEner trna+ f = rnaFoldMFEEner frna+{-# Inlinable mfeDelta' #-}++centroidDelta' :: Bool → Bool → BM.ScaleFunction Double+centroidDelta' mx sq frna trna = guardPositives mx sq $ t - f+ where t = rnaFoldCentroidEner trna+ f = rnaFoldCentroidEner frna+{-# Inlinable centroidDelta' #-}++guardPositives mx sq = (if sq then (\z → if z > 0 then z^2 else z) else id) . (if mx then max 0 else id)++mfebpdist' ∷ Bool → Bool → BM.ScaleFunction Double+mfebpdist' mx sq frna trna = fromIntegral $ d1Distance (bldD1S $ rnaFoldMFEStruc frna) (bldD1S $ rnaFoldMFEStruc trna)++centroidbpdist' ∷ Bool → Bool → BM.ScaleFunction Double+centroidbpdist' mx sq frna trna = fromIntegral $ d1Distance (bldD1S $ rnaFoldCentroidStruc frna) (bldD1S $ rnaFoldCentroidStruc trna)++partfun' ∷ BM.ScaleFunction Double → BM.ScaleFunction (Log Double)+partfun' f frna trna = Exp . negate $ f frna trna+{-# Inlineable partfun' #-}++rna2dna ∷ RNAfoldResult → RNAfoldResult+rna2dna r = r { rnaFoldSequence = BS.map go $ rnaFoldSequence r } where+ go x =let x' = toUpper x in if x' == 'U' then 'T' else x'++-- | Given ancestral and extant sequence, and possibly an intermediate+-- mutation, as well as a list of intermediates create the set to RNA structure+-- mapping.+--+-- TODO some of the functions here should go into @Lib-SequencePolymorphism@.+--+-- TODO should run within @ExceptT@ !++genSet+ ∷ Ancestral+ → Extant+ → Maybe (Int,Char)+ -- ^ If @Just@ then the 0-based position and character of the intermediate+ -- mutation.+ → [RNAfoldResult]+ → (HM.HashMap Int RNAfoldResult, Int, Int, B.BimapHashMap Int Int)+genSet (Ancestral a') (Extant e') v xs = (HM.fromList kv, ipos, B.size posbit, posbit)+ where kv = [ (b, maybe (kvErr b) id $ HM.lookup (pat2str b) lkupRes) | b ← bits ]+ kvErr b = error $ show (lkupRes, sort $ map rnaFoldSequence xs, v, posbit, b, bits, pat2str b)+ -- update a/e based on if we have the intermediate mutation set up.+ a = maybe a' (\(i,c) → unpackedChars.ix i .~ c $ a') v+ e = maybe e' (\(i,c) → unpackedChars.ix i .~ c $ e') v+ -- all positions where the two bytestrings differ, together with the+ -- differing characters+ ks = filter (\(_,i,j) → i/=j) $ zip3 [0∷Int ..] (BS.unpack a') (BS.unpack e')+ -- turn into a bijection of actual position (first) and bit in bitset+ -- (second)+ posbit ∷ B.BimapHashMap Int Int+ posbit = B.fromList $ zip (ks^..traverse._1) [0∷Int ..]+ -- all bit patterns+ bits = [0 .. 2^B.size posbit - 1]+ -- convert a bit pattern to an actual string, to be looked up. Start+ -- with the ancestral sequence and for each @1@, modify the character+ -- into the one encountered in the extant sequence.+ pat2str ∷ Int → ByteString+ pat2str = let go s k = unpackedChars.ix (lk k) .~ (e `BS.index` lk k) $ s+ lk = maybe (error "lk") id . B.lookupR posbit+ in foldl' go a . activeBitsL+ -- lookup from sequence to RNAfoldResult+ lkupRes = HM.fromList [ (rnaFoldSequence x,x) | x ← xs ]+ -- is the global mutation intermediate? Yes: then >= 0 is the bitset+ -- element to be returned here+ ipos = maybe (-1) (\(z,_) → if BS.index a' z /= BS.index e' z then (maybe (error "ipos") id $ B.lookupL posbit z) else -1) v++-- | Run the intermediate / backmutation order variant. This variant is slow,+-- and requires large pre-calculated files, we parallelize and aggregate as+-- much as possible.+--+-- TODO read monad ?!++{-+runBackmutationVariants+ ∷ Int+ → [Char]+ → GlobalBackmutations+ → [BackmutationCol]+ → Ancestral+ → Extant+ → ExceptT String IO ()+runBackmutationVariants aggregate alphabet globback backcols ancestral extant = do+ -- Load all sequences for the original problem -- they are needed anyway+ (seqCount, origSeqs', variants) ← createRNAlandscape2 alphabet globback backcols ancestral extant+ let origSeqs = Trie.fromList [ (s,()) | s ← origSeqs' ]+ origStrs ← filter (\r → rnaFoldSequence r `Trie.member` origSeqs) <$> readRNAfoldFiles (error "workdb")+ let rnas = error "bitset -> rnafoldresult data"+ -- Group into sets of @aggregate@ elements for sequence aggregation+ let ass = chunksOf aggregate variants+ forM_ ass $ \as → do+ -- the required sequences are given by @origSeqs@ but modified at the appropriate position+ let aggrSeqss = map (\(p,n,xs) → (p,n,Trie.fromList [ (x,()) | x ← xs ])) as+ let allss = foldl' (\z (_,_,x) → Trie.unionL z x) Trie.empty aggrSeqss+ -- read in the structures for all as+ aggrStr ← filter (\r → rnaFoldSequence r `Trie.member` allss) <$> readRNAfoldFiles (error "workdb")+ let go (p,n,xs) = (p,n) where+ -- prepare the @ntrs@ data structure for each as+ ntrs = undefined $ filter (\r → rnaFoldSequence r `Trie.member` xs) aggrStr+ fMnD = BM.forwardMinDist (error "number of known mutations") (error "scale function") rnas ntrs+ fwdZ = BM.forwardEvidence (error "number of known mutations") (error "scale function for evidence") rnas ntrs+ -- parallel calculation and output for the different cases+ let rs = Par.parMap Par.rdeepseq go aggrSeqss+ forM_ rs $ \r → do+ return ()+ return ()+ -- print each output+ return ()+-}
+ BioInf/MutationOrder/BackMutations.hs view
@@ -0,0 +1,390 @@++-- | Specialized mutation order grammar and algebras that incorporate exactly+-- one intermediate mutation. Not necessarily only of the backmutation kind.+--+-- Calculated are (i) the total evidence @Z@ for two input sequences. (ii) The+-- minimal weight distance.++module BioInf.MutationOrder.BackMutations where++import qualified Data.Char as C+import Data.Bits+import Data.HashMap.Strict (HashMap)+import qualified Data.HashMap.Strict as HM+import qualified Data.Vector.Fusion.Stream.Monadic as SM+import Numeric.Log+import Data.Text (Text)+import qualified Data.Text as T+import qualified Data.ByteString.Char8 as BS+import Text.Printf+import qualified Data.Vector.Unboxed as VU++import qualified Data.Bijection.HashMap as B+import Data.PrimitiveArray hiding (toList,map)+import qualified Data.PrimitiveArray as PA+import FormalLanguage+import ADP.Fusion.Core hiding (PeekIndex)+import ADP.Fusion.Set1 hiding (PeekIndex)+import ADP.Fusion.Term.PeekIndex.Set1+import ADP.Fusion.Unit hiding (PeekIndex)++import BioInf.MutationOrder.SequenceDB (RNAfoldResult(..))++[formalLanguage|+Verbose+Grammar: MinDist+N: Begin+N: Inter+N: Final+N: Start+T: peek -- insert backmutation+T: edg -- normal edge+T: nd -- single node+S: Start+-- A number of mutations that do not include the intermediate mutation.+Begin -> beginEmpty <<< ε+Begin -> beginNode <<< nd+Begin -> beginEdge <<< Begin edg+-- Insert the single intermediate mutation, followed by more known events.+-- TODO this *could* be done with a statically active guard!+Inter -> interFromBegin <<< Begin peek -- activate backmutation or intermediate mutation+Inter -> interEdge <<< Inter edg+-- Undo the intermediate mutation, followed by more known events.+Final -> finalFromInter <<< Inter edg -- this moves from an intermediate to a final mutation+Final -> finalUnBackmut <<< Inter peek -- deactivates the backmutation+Final -> finalEdge <<< Final edg+Start -> finis <<< Final+//+Emit: MinDist+|]++makeAlgebraProduct ''SigMinDist++++type ScaleFunction d = RNAfoldResult → RNAfoldResult → d++-- | Minimal distance calculation under the influence of one intermediate+-- mutation.+--+-- NOTE be very careful to check that @rnas@ (no intermediate mutation) and+-- @ntrs@ are used correctly!++aMinDist+ ∷ Monad m+ ⇒ d+ -- ^ neutral element for @omin@.+ → (d → d → d)+ -- ^ omin+ → (d → d → d)+ -- ^ oplus+ → ScaleFunction d+ -- ^ Combine two 'RNAfoldResult's into the resulting @d@ score.+ → Int+ -- ^ Index position of the intermediate mutation or @(-1)@ if independent of the observed mutations.+ -- TODO IN COORDINATES OF THE BITSET SPACE!+ → HashMap Int RNAfoldResult+ -- ^ RNAs without the intermediate mutation set+ → HashMap Int RNAfoldResult+ -- ^ RNAs with the intermediate mutation set+ → SigMinDist m d d (Int:.From:.To) (Int:.To) (BS1 Last I)+aMinDist neutral omin oplus scaled ipos rnas ntrs = SigMinDist+ { beginEmpty = \() → neutral+ -- ^ Not a single mutation has happened.+ , beginNode = \(nset:.To n) → let frna = rnas HM.! 0+ trna = rnas HM.! (0 `setBit` n)+ in if ipos < 0 || n /= ipos+ then scaled frna trna+ else neutral+ -- ^ Activate a single, first mutation.+ , beginEdge = \x (fset:.From f:.To t) → let frna = rnas HM.! fset+ trna = rnas HM.! tset+ tset = fset `setBit` f `setBit` t+ in if ipos < 0 || (not $ tset `testBit` ipos)+ then x `oplus` scaled frna trna+ else neutral+ -- ^ These edges are mutational events without the influence of the+ -- intermediate mutation present.+ , interFromBegin = \x (BS1 fset b) → let frna = rnas HM.! getBitSet fset+ trna = ntrs HM.! getBitSet fset+ -- Only allow if either the mutation does not+ -- influence the observed mutations or the+ -- influenced mutation is still available.+ in if ipos < 0 || (not $ fset `testBit` ipos)+ then x `oplus` scaled frna trna+ else neutral+ -- ^ The intermediate mutation is about to be set. Only if @ipos@ is+ -- independent of the other mutational events or we have not tried activating it.+ , interEdge = \x (fset:.From f:.To t) → let frna = ntrs HM.! fset+ trna = ntrs HM.! tset+ tset = fset `setBit` f `setBit` t+ in if ipos < 0 || (not $ tset `testBit` ipos)+ then x `oplus` scaled frna trna+ else neutral+ -- ^ These are edges inserted while the intermediate mutation is active. Only+ -- allow intermediate mutations that do not move the @ipos@ mutation into its+ -- final state.+ -- TODO check if @f@ or @t@ will be set here+ , finalFromInter = \x (fset:.From f:.To t) → let frna = ntrs HM.! fset+ trna = rnas HM.! tset+ tset = fset `setBit` f `setBit` t+ in if ipos == t+ then x `oplus` scaled frna trna+ else neutral+ -- ^ Flips the intermediate mutation (@ipos >= 0@) into the final state.+ -- TODO check if @t@ is the one being set+ , finalUnBackmut = \x (BS1 fset t) → let frna = ntrs HM.! getBitSet fset+ trna = rnas HM.! getBitSet fset+ -- Only allow if either the mutation does not+ -- influence the observed mutations or the+ -- influenced mutation is + in if ipos < 0+ then x `oplus` scaled frna trna+ else neutral+ -- ^ Now the intermediate mutation is undone. Is only used if @ipos@ is+ -- independent (@==(-1)@) of the observed events.+ , finalEdge = \x (fset:.From f:.To t) → let frna = rnas HM.! fset+ trna = rnas HM.! tset+ tset = fset `setBit` f `setBit` t+ in if ipos < 0 || (fset `testBit` ipos)+ then x `oplus` scaled frna trna+ else neutral+ , finis = id+ -- ^ Collapse over possible end points+ , h = SM.foldl' omin neutral+ -- ^ Find the minimal distance from ancestral to extant sequence.+ }+{-# Inline aMinDist #-}++-- |+--+-- TODO Use text builder++aPretty+ ∷ Monad m+ ⇒ ScaleFunction Double+ → Int+ -- ^ index of intermediate mutation+ → Int+ -- ^ real ipos value+ → HashMap Int RNAfoldResult+ -- ^ rnas+ → HashMap Int RNAfoldResult+ -- ^ intermediate rnas+ → B.BimapHashMap Int Int+ -- ^ actual mutation position / bit in bitset+ → SigMinDist m Text [Text] (Int:.From:.To) (Int:.To) (BS1 Last I)+aPretty scaled ipos realPos rnas ntrs mutpos = SigMinDist+ { beginEmpty = \() → ""+ , beginNode = \(nset:.To n) → let frna = rnas HM.! 0+ trna = rnas HM.! (0 `setBit` n)+ in if ipos < 0 || n /= ipos+ then T.concat [showHdr "Ancestral" frna, showMut "1st Mutation" frna trna n]+ else ""+ -- ^ Activate a single, first mutation.+ , beginEdge = \x (fset:.From f:.To t) → let frna = rnas HM.! fset+ trna = rnas HM.! tset+ tset = fset `setBit` f `setBit` t+ in if ipos < 0 || ipos /= t+ then T.concat [x, showMut "Edge (Pre)" frna trna t]+ else ""+ -- ^ These edges are mutational events without the influence of the+ -- intermediate mutation present.+ , interFromBegin = \x (BS1 fset b) → let frna = rnas HM.! getBitSet fset+ trna = ntrs HM.! getBitSet fset+ imorbm = if ipos < 0 then "Activate BM" else "Activate IM"+ -- Only allow if either the mutation does not+ -- influence the observed mutations or the+ -- influenced mutation is still available.+ in if ipos < 0 || (not $ fset `testBit` ipos)+ then T.concat [x, showBackmut imorbm frna trna realPos]+ else ""+ -- ^ The intermediate mutation is about to be set. Only if @ipos@ is+ -- independent of the other mutational events or we have not tried activating it.+ , interEdge = \x (fset:.From f:.To t) → let frna = ntrs HM.! fset+ trna = ntrs HM.! tset+ tset = fset `setBit` f `setBit` t+ in if ipos < 0 || (not $ tset `testBit` ipos)+ then T.concat [x, showMut "Edge (Active)" frna trna t]+ else ""+ -- ^ These are edges inserted while the intermediate mutation is active. Only+ -- allow intermediate mutations that do not move the @ipos@ mutation into its+ -- final state.+ -- TODO check if @f@ or @t@ will be set here+ , finalFromInter = \x (fset:.From f:.To t) → let frna = ntrs HM.! fset+ trna = rnas HM.! tset+ tset = fset `setBit` f `setBit` t+ in if ipos == t+ then T.concat [x, showUnBackmut "Deactivate IM" frna trna realPos]+ else ""+ -- ^ Flips the intermediate mutation (@ipos >= 0@) into the final state.+ -- TODO check if @t@ is the one being set+ , finalUnBackmut = \x (BS1 fset t) → let frna = ntrs HM.! getBitSet fset+ trna = rnas HM.! getBitSet fset+ -- Only allow if either the mutation does not+ -- influence the observed mutations or the+ -- influenced mutation is + in if ipos < 0+ then T.concat [x, showUnBackmut "Deactivate BM" frna trna realPos]+ else ""+ -- ^ Now the intermediate mutation is undone. Is only used if @ipos@ is+ -- independent (@==(-1)@) of the observed events.+ , finalEdge = \x (fset:.From f:.To t) → let frna = rnas HM.! fset+ trna = rnas HM.! tset+ tset = fset `setBit` f `setBit` t+ in if ipos < 0 || (fset `testBit` ipos)+ then T.concat [x, showMut "Edge (Post)" frna trna t]+ else ""+ , finis = id+ -- ^ Collapse over possible end points+ , h = SM.toList . SM.filter (not . T.null)+ -- ^ Find the minimal distance from ancestral to extant sequence.+ } where+ muts = let as = VU.generate n (const ' ')+ n = BS.length $ rnaFoldSequence $ rnas HM.! 0+ bs = as VU.// [ (k*10 - 1, C.chr $ C.ord '0' + k `mod` 10) | k ← [1 .. n `div` 10] ]+ cs = bs VU.// [ (k, 'v') | (k,_) ← B.toList mutpos ]+ ds = cs VU.// [ if ipos < 0 then (realPos-1,'!') else (realPos-1,'+') ]+ in VU.toList ds -- ++ show (ipos,realPos)+ showHdr ∷ String → RNAfoldResult → Text+ showHdr what frna = T.pack $ printf "%s\n%s%s %s\n%s%s %5.1f\n%s%s %5.1f\n"+ (spaces ++ muts)+ spaces+ (BS.unpack $ rnaFoldSequence frna)+ what+ spaces+ (BS.unpack $ rnaFoldMFEStruc frna)+ (rnaFoldMFEEner frna)+ spaces+ (BS.unpack $ rnaFoldCentroidStruc frna)+ (rnaFoldCentroidEner frna)+ showMut ∷ String → RNAfoldResult → RNAfoldResult → Int → Text+ showMut what frna trna p = T.pack $ printf "%5d %5.1f %s %s\n%s%s %5.1f\n%s%s %5.1f\n"+ (maybe (error $ "showMut" ++ show p) (+1) $ B.lookupR mutpos p)+ (deltaE frna trna)+ (BS.unpack $ rnaFoldSequence trna)+ what+ spaces+ (BS.unpack $ rnaFoldMFEStruc trna)+ (rnaFoldMFEEner trna)+ spaces+ (BS.unpack $ rnaFoldCentroidStruc trna)+ (rnaFoldCentroidEner trna)+ showBackmut ∷ String → RNAfoldResult → RNAfoldResult → Int → Text+ showBackmut what frna trna p = T.pack $ printf "%c %3d %5.1f %s %s\n%s%s %5.1f\n%s%s %5.1f\n"+ (if ipos < 0 then '!' else '+')+ p+ (deltaE frna trna)+ (BS.unpack $ rnaFoldSequence trna)+ what+ spaces+ (BS.unpack $ rnaFoldMFEStruc trna)+ (rnaFoldMFEEner trna)+ spaces+ (BS.unpack $ rnaFoldCentroidStruc trna)+ (rnaFoldCentroidEner trna)+ showUnBackmut ∷ String → RNAfoldResult → RNAfoldResult → Int → Text+ showUnBackmut what frna trna p = T.pack $ printf "%c %3d %5.1f %s %s\n%s%s %5.1f\n%s%s %5.1f\n"+ (if ipos < 0 then '!' else '+')+ p+ (deltaE frna trna)+ (BS.unpack $ rnaFoldSequence trna)+ what+ spaces+ (BS.unpack $ rnaFoldMFEStruc trna)+ (rnaFoldMFEEner trna)+ spaces+ (BS.unpack $ rnaFoldCentroidStruc trna)+ (rnaFoldCentroidEner trna)+ deltaE ∷ ScaleFunction Double+ deltaE = scaled -- frna trna = rnaFoldMFEEner trna - rnaFoldMFEEner frna+ spaces = replicate 12 ' '+{-# Inline aPretty #-}++type FwdBS1 x = TwITbl Id Unboxed EmptyOk (BS1 Last I) x+type FwdUnit x = TwITbl Id Unboxed EmptyOk (Unit I) x++type BTS1 x b = TwITblBt Unboxed EmptyOk (BS1 Last I) x Id Id b+type BTUnit x b = TwITblBt Unboxed EmptyOk (Unit I) x Id Id b++-- |+--+-- TODO check if inlining of @ScaleFunction@ improves performance+-- substantially.++forwardMinDist+ ∷ Int+ → ScaleFunction Double+ → Int+ → HashMap Int RNAfoldResult+ → HashMap Int RNAfoldResult+ → Z:.FwdBS1 Double:.FwdBS1 Double:.FwdBS1 Double:.FwdUnit Double+forwardMinDist n scaled ipos rnas ntrs =+ let + in mutateTablesST $ gMinDist (aMinDist 999999 min (+) scaled ipos rnas ntrs)+ (ITbl 0 0 EmptyOk (fromAssocs (BS1 0 (-1)) (BS1 (2^n-1) (Boundary $ n-1)) 999999 [])) -- Begin+ (ITbl 2 0 EmptyOk (fromAssocs (BS1 0 (-1)) (BS1 (2^n-1) (Boundary $ n-1)) 999999 [])) -- Final+ (ITbl 1 0 EmptyOk (fromAssocs (BS1 0 (-1)) (BS1 (2^n-1) (Boundary $ n-1)) 999999 [])) -- Inter+ (ITbl 3 0 EmptyOk (fromAssocs Unit Unit 999999 [])) -- Start+ EdgeWithSet -- normal mutational event+ Singleton -- first mutational event+ (PeekIndex ∷ PeekIndex (BS1 Last I)) -- undo intermediate+{-# NoInline forwardMinDist #-}++forwardMinDistValue fwd = md PA.! PA.Unit+ where (Z:.fwdB:.fwdF:.fwdI:.fwdS) = fwd+ TW (ITbl _ _ _ md) _ = fwdS++backtrackMinDist1+ ∷ Int+ → ScaleFunction Double+ → Int+ → Int+ → HashMap Int RNAfoldResult+ → HashMap Int RNAfoldResult+ → B.BimapHashMap Int Int+ → Z:.FwdBS1 Double:.FwdBS1 Double:.FwdBS1 Double:.FwdUnit Double+ → [Text]+backtrackMinDist1 n scaled ipos realPos rna ntrs mutpos (Z:.fwdB:.fwdF:.fwdI:.fwdS) = unId $ axiom btS+ where !(Z:.btB:.btF:.btI:.btS) =+ gMinDist (aMinDist 999999 min (+) scaled ipos rna ntrs <|| aPretty scaled ipos realPos rna ntrs mutpos)+ (toBacktrack fwdB (undefined :: Id a -> Id a))+ (toBacktrack fwdF (undefined :: Id a -> Id a))+ (toBacktrack fwdI (undefined :: Id a -> Id a))+ (toBacktrack fwdS (undefined :: Id a -> Id a))+ EdgeWithSet+ Singleton+ (PeekIndex ∷ PeekIndex (BS1 Last I))+ :: Z:.BTS1 Double Text:.BTS1 Double Text:.BTS1 Double Text:.BTUnit Double Text+{-# NoInline backtrackMinDist1 #-}++-- |+--+-- TODO check if inlining of @ScaleFunction@ improves performance+-- substantially.++forwardEvidence+ ∷ Int+ → ScaleFunction (Log Double)+ → Int+ → HashMap Int RNAfoldResult+ → HashMap Int RNAfoldResult+ → Z:.FwdBS1 (Log Double):.FwdBS1 (Log Double):.FwdBS1 (Log Double):.FwdUnit (Log Double)+forwardEvidence n scaled ipos rnas ntrs =+ let + in mutateTablesST $ gMinDist (aMinDist 0 (+) (*) scaled ipos rnas ntrs)+ (ITbl 0 0 EmptyOk (fromAssocs (BS1 0 (-1)) (BS1 (2^n-1) (Boundary $ n-1)) 0 [])) -- Begin+ (ITbl 2 0 EmptyOk (fromAssocs (BS1 0 (-1)) (BS1 (2^n-1) (Boundary $ n-1)) 0 [])) -- Final+ (ITbl 1 0 EmptyOk (fromAssocs (BS1 0 (-1)) (BS1 (2^n-1) (Boundary $ n-1)) 0 [])) -- Inter+ (ITbl 3 0 EmptyOk (fromAssocs Unit Unit 0 [])) -- Start+ EdgeWithSet -- normal mutational event+ Singleton -- first mutational event+ (PeekIndex ∷ PeekIndex (BS1 Last I)) -- undo intermediate+{-# NoInline forwardEvidence #-}++forwardEvidenceValue fwd = md PA.! PA.Unit+ where (Z:.fwdB:.fwdF:.fwdI:.fwdS) = fwd+ TW (ITbl _ _ _ md) _ = fwdS+++
BioInf/MutationOrder/MinDist.hs view
@@ -121,7 +121,7 @@ -- -- TODO Use text builder -aPretty :: Monad m => ScaleFunction -> Landscape -> SigMinDist m Text [Text] (Int:.From:.To) (Int:.To)+aPretty :: Monad m => ScaleFunction -> Landscape -> SigMinDist m (Text,[Int]) [(Text,[Int])] (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)@@ -130,8 +130,8 @@ 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 = \() -> ""+ in (T.concat [fst x, showMut frna trna t' eM eC eS], t' : snd x)+ , mpty = \() -> ("",[]) , node = \(nset:.To n) -> let frna = rnas HM.! (BitSet 0)@@ -140,7 +140,7 @@ 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]+ in (T.concat [showHdr frna n', showMut frna trna n' eM eC eS], [n']) , fini = id , h = SM.toList } where@@ -175,14 +175,14 @@ -type TS1 x = TwITbl Id Unboxed EmptyOk (BS1 First I) x+--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 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 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@@ -204,14 +204,14 @@ Singleton {-# NoInline forwardMinDist1 #-} -backtrackMinDist1 :: ScaleFunction -> Landscape -> Z:.TS1L Double:.U Double -> [Text]+backtrackMinDist1 :: ScaleFunction -> Landscape -> Z:.TS1L Double:.U Double -> [(Text,[Int])] 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+ :: Z:.BT1L Double (Text,[Int]):.BTU Double (Text,[Int]) {-# NoInline backtrackMinDist1 #-} -- | Count the number of co-optimals@@ -242,7 +242,7 @@ -- 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 -> (Double,[(Text,[Int])]) runCoOptDist scaleFunction landscape = (unId $ axiom fwdu,bs) where !(Z:.fwd1:.fwdu) = forwardMinDist1 scaleFunction landscape bs = backtrackMinDist1 scaleFunction landscape (Z:.fwd1:.fwdu)@@ -255,21 +255,21 @@ -- | 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 #-}+--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 =
BioInf/MutationOrder/RNA.hs view
@@ -10,30 +10,35 @@ module BioInf.MutationOrder.RNA where -import Data.Aeson as DA-import Data.Bits+import Data.Char (toUpper) import Codec.Compression.GZip (compress,decompress) import Control.Arrow (second) import Control.DeepSeq+import Control.Error+import Control.Monad (unless) import Control.Parallel.Strategies+import Data.Aeson as DA+import Data.Bits import Data.ByteString (ByteString)+import Data.Char (isDigit)+import Data.List (sort,nub,(\\)) import Data.Maybe (catMaybes)+import Data.Monoid import Data.Serialize+import Data.Serialize.Instances+import Data.Text.Encoding (decodeUtf8, encodeUtf8)+import Data.Tuple (swap) 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.HashMap.Strict as HM+import qualified Data.IntMap.Strict as IM 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@@ -193,6 +198,77 @@ landscapeDestination <- encodeUtf8 <$> v .: "landscapeDestination" mutationPositions <- v .: "mutationPositions" return Landscape{..}++newtype Ancestral = Ancestral { getAncestral ∷ ByteString }++newtype Extant = Extant { getExtant ∷ ByteString }++newtype GlobalBackmutations = GlobalBackmutations { getGlobalBackmutations ∷ Int }++newtype BackmutationCol = BackmutationCol Int++-- | The new method for creating the RNA landscape. The outgoing stream of+-- sequences should be directly to the 'BioInf.MutationOrder.SequenceDB'+-- module, and 'writeSequenceFiles' in particular.+--+-- TODO generalise over the alphabet+--+-- TODO currently, the number of backmutations is either @0@ or @1@. For more,+-- a number of things need to be fixed: (i) creation of all variants of+-- @1,2,3...g@ backmutations. That is, consider all sequences that have at+-- least one, but no more than @g@ global back mutations.+--+-- TODO this will create duplicate sequences for those columns where the+-- ancestral and extant have different nucleotides during global backmutation+-- generation.++createRNAlandscape2+ ∷ (Monad m)+ ⇒ [Char]+ → Either GlobalBackmutations [Int]+ → [BackmutationCol]+ → Ancestral+ → Extant+ → ExceptT String m (Int,[ByteString], [(Int,Char,[ByteString])])+ -- ^ Return A complex triple with (total sequence count, original problem+ -- sequences, list of global variants). The list of global variants holds for+ -- each (index,single nucleotide intermediate, list of variant sequences).+createRNAlandscape2 alphabet gxorgs bs (Ancestral a) (Extant e) = do+ -- some sanity checks+ unless (BS.length a == BS.length e) $ throwE "different sequence lengths encountered"+ -- expand later! Right now, we either generate the + unless (either ((<=1) . getGlobalBackmutations) ((<=1) . length) gxorgs)+ $ throwE "we currently allow *at most* one globally active backmutation"+ -- TODO need to fix up interesting columns+ unless (null bs) $ throwE "fix up interesting columns"+ -- collect the possible characters for each position.+ let ahm = IM.fromList . zip [0∷Int ..] . map ((:[]) . toUpper) $ BS.unpack a+ let ehm = IM.fromList . zip [0..] . map ( (:[]) . toUpper) $ BS.unpack e+ -- back mutation columns are active together with the above+ let bhm = IM.fromListWith (++) [ (k,alphabet) | BackmutationCol k ← bs ]+ let merge x y = sort . nub $ x++y+ let hm = IM.unionWith merge (IM.unionWith merge ahm ehm) bhm+ -- global back mutations, this will introduce only those characters not already present in each column+ let gbm = IM.fromList [ (k, alphabet \\ hm IM.! k)+ | k ← either (\_ → [0..BS.length a -1]) id gxorgs ]+ -- begin with the set of sequences without any global backmutations+ let localList = map BS.pack . sequence . map snd . IM.toAscList+ let localCount = product . map (length . snd) $ IM.toAscList hm -- do *not* count explicitly!+ let globalModifiers = concat . map (\(k,xs) → map (k,) xs) $ IM.toAscList gbm+ -- we now repeat the local generation, but "splice in" the globally modified characters+ -- TODO we currently allow exactly one global modifier in @[(k,z)]@+ let globals = if True -- g == 1 -- TODO need to fix up+ then [ (k, z, [ BS.take k orig `BS.append` (BS.cons z (BS.drop (k+1) orig))+ | orig ← localList hm+ -- produce a sequence only if this mutation is+ -- independent of the a/e mutations OR we look+ -- at a variant or the ancestral sequence. This+ -- prevents duplicate entries.+ , (a `BS.index` k == e `BS.index` k) || (a `BS.index` k == orig `BS.index` k) ] )+ | (k,z) ← globalModifiers ]+ else []+ let globalCount = localCount * length globalModifiers -- if g == 1 then localCount * length globalModifiers else 0+ return $ (localCount + globalCount, localList hm, globals) -- | --
+ BioInf/MutationOrder/SequenceDB.hs view
@@ -0,0 +1,167 @@++-- | Functions to deal with a (large) RNA sequence data base and corresponding+-- RNAfold foldings.++module BioInf.MutationOrder.SequenceDB where++import Codec.Compression.GZip (compress,decompress)+import Control.Error+import Control.Monad.Except+import Control.Monad.IO.Class (liftIO, MonadIO)+import Data.ByteString.Char8 (ByteString)+import Data.Char (isDigit)+import Data.Hashable+import Debug.Trace+import GHC.Generics+import qualified Data.Attoparsec.ByteString.Char8 as AC+import qualified Data.Attoparsec.ByteString.Lazy as A+import qualified Data.ByteString.Char8 as BS+import qualified Data.ByteString.Lazy.Char8 as BSL+import qualified Data.Char as C+import System.Directory (doesFileExist, getDirectoryContents)+import System.FilePath ((</>), (<.>))+import System.IO.Unsafe (unsafeInterleaveIO)+import qualified System.FilePath.Find as FP++++type PrefixLen = Int+type SeqsPerFile = Int++-- * Plain sequences, with no folding yet++-- | Writes sequence files out to disk. Does *not* check if duplicates are+-- written. The user should do this externally via @sort@ and @uniq@. The next+-- free prefix is chosen.++writeSequenceFiles+ ∷ (MonadIO m)+ ⇒ FilePath+ -- ^ Directory to write to.+ → PrefixLen+ -- ^ Number of characters in the prefix to count up. Fail if the prefix space+ -- is too small.+ → SeqsPerFile+ -- ^ How many sequences to write per file+ → [ByteString]+ -- ^ Sequences to write out.+ → ExceptT String m ()+writeSequenceFiles fp pfx spf xs = do+ -- cull prefix list to unused prefixes only+ let unused [] = throwE "writeSequenceFiles: prefix space empty"+ unused (x:xs) = do+ dfe ← liftIO . doesFileExist $ fp </> x <.> "gz"+ if dfe+ then unused xs+ else return $ x:xs+ -- will write out file names to unused prefixes.+ let go ps [] = return ()+ go ps ys = do+ let (here,there) = splitAt spf ys+ (p:qs) ← unused ps+ -- write to compressed file @p@+ let fname = fp </> p <.> "gz"+ liftIO . BSL.writeFile fname . compress . BSL.unlines $ map BSL.fromStrict here+ go qs there+ let pfxs = sequence $ replicate pfx ['a'..'z']+ go pfxs xs++-- | Reads all sequences.+--+-- TODO make sure to read lazily enough that collection happens without first+-- reading in the whole list.++readSequenceFiles+ ∷ FilePath+ → ExceptT () IO [ByteString]+readSequenceFiles fp = do+ return undefined++-- * Structural information from RNAfold++-- | Space-efficient RNAfold structure++data RNAfoldResult = RNAfoldResult+ { rnaFoldSequence ∷ !ByteString+ , rnaFoldMFEStruc ∷ !ByteString+ , rnaFoldMFEEner ∷ !Double+ , rnaFoldMfeFrequency ∷ !Double+ , rnaFoldEnsembleStruc ∷ !ByteString+ -- ^ uses special syntax with unpaired, weakly paired, somewhat paired,+ -- somewhat paired up or down, strongly paired up or down for the ensemble+ , rnaFoldEnsembleEner ∷ !Double+ -- ^ this *is* the ensemble free energy+ , rnaFoldCentroidStruc ∷ !ByteString+ , rnaFoldCentroidEner ∷ !Double+ , rnaFoldDiversity ∷ !Double+ }+ deriving (Read,Show,Eq,Ord,Generic)++-- | Lazily read @RNAfold@ structures.+--+-- TODO use @pipes/machines@! we need lazy reading of files and live in an exceptt transformer stack!+-- TODO generalize transformer stack++readRNAfoldFiles+ ∷ FilePath+ → IO [RNAfoldResult]+readRNAfoldFiles fp = do+ let go [] = return []+ go (f:fs) = do+ bs ← unsafeInterleaveIO (putStrLn ("#" ++ f) >> decompress <$> BSL.readFile f)+ let rs = either error id $ runExcept (bslToRNAfoldResult bs)+ rss ← go fs+ return $ rs ++ rss+ go ∷ [FilePath] → IO [RNAfoldResult]+ FP.find FP.always (FP.extension FP.==? ".gz") (fp </> "structures") >>= go+{-# NoInline readRNAfoldFiles #-}++-- |++bslToRNAfoldResult ∷ (Monad m) ⇒ BSL.ByteString → ExceptT String m [RNAfoldResult]+bslToRNAfoldResult bs = do+ case A.eitherResult $ A.parse pRNAfold bs of+ Left e → throwE e+ Right r → return r+{-# Inline bslToRNAfoldResult #-}++-- |+--+-- @+-- echo "CCCAAAGGG\nCCCAAAGGG" | ./RNAfold -p+-- CCCAAAGGG+-- (((...))) ( -1.20)+-- (((...))) [ -1.41]+-- (((...))) { -1.20 d=1.06}+-- frequency of mfe structure in ensemble 0.707288; ensemble diversity 1.67 +-- CCCAAAGGG+-- (((...))) ( -1.20)+-- (((...))) [ -1.41]+-- (((...))) { -1.20 d=1.06}+-- frequency of mfe structure in ensemble 0.707288; ensemble diversity 1.67 +-- @++pRNAfold ∷ A.Parser [RNAfoldResult]+pRNAfold = A.many1' go <* A.endOfInput where+ go = do+ -- 1. sequence+ rnaFoldSequence ← BS.copy <$> AC.takeWhile AC.isAlpha_ascii <* AC.skipSpace A.<?> "RNAfold sequence"+ -- 2. mfe+ rnaFoldMFEStruc ← BS.copy <$> AC.takeTill AC.isSpace <* AC.skipSpace A.<?> "RNAfold MFE structure"+ rnaFoldMFEEner ← AC.char '(' *> AC.skipSpace *> AC.double <* AC.char ')' <* AC.skipSpace A.<?> "RNAfold MFE energy"+ -- 3. ensemble+ rnaFoldEnsembleStruc ← BS.copy <$> AC.takeTill AC.isSpace <* AC.skipSpace+ rnaFoldEnsembleEner ← AC.char '[' *> AC.skipSpace *> AC.double <* AC.char ']' <* AC.skipSpace+ -- 4. centroid+ rnaFoldCentroidStruc ← BS.copy <$> AC.takeTill AC.isSpace <* AC.skipSpace+ rnaFoldCentroidEner ← AC.char '{' *> AC.skipSpace *> AC.double <* AC.skipSpace+ dequal ← AC.string "d=" *> AC.double <* AC.char '}' <* AC.skipSpace+ -- 5.mfe frequency and diversity+ AC.string "frequency of mfe structure in ensemble" *> AC.skipSpace A.<?> "frequency"+ rnaFoldMfeFrequency ← AC.double+ AC.string "; ensemble diversity" *> AC.skipSpace+ rnaFoldDiversity ← AC.double+ AC.skipSpace+ return RNAfoldResult{..}+{-# Inline pRNAfold #-}+
MutationOrder.cabal view
@@ -1,7 +1,7 @@ name: MutationOrder-version: 0.0.0.2+version: 0.0.1.0 author: Maria Beatriz Walter Costa, Christian Hoener zu Siederdissen, 2017-copyright: Christian Hoener zu Siederdissen, 2017+copyright: Maria Beatriz Walter Costa, 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@@ -11,19 +11,23 @@ build-type: Simple stability: experimental cabal-version: >= 1.10.0-tested-with: GHC == 7.10.3, GHC == 8.0.1+tested-with: GHC == 8.0.2 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.+ Developed to analyse the @HAR 1@ region, but agnostic to the+ actual sequences and can be used to analyze any RNA sequence+ that fits the algorithmic constraints. . 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.+ .+ We currently only consider point mutations, not in-dels. @@ -48,16 +52,24 @@ library build-depends: base >= 4.7 && < 5.0 , aeson >= 1.1+ , attoparsec >= 0.13 , bytestring+ , bytestring-trie >= 0.2 , cereal >= 0.5 , cereal-vector >= 0.2 , containers , deepseq >= 1.4 , directory+ , errors >= 2.0+ , filemanip >= 0.3 , filepath+ , hashable >= 1.2+ , lens >= 4.0 , log-domain >= 0.10+ , mtl , parallel >= 3.2 , serialize-instances >= 0.1+ , split >= 0.2 , text >= 1.0 , unordered-containers >= 0.2.7 , vector >= 0.11@@ -70,6 +82,7 @@ , BiobaseXNA == 0.9.3.* , DPutils == 0.0.1.* , FormalGrammars == 0.3.1.*+ , OrderedBits == 0.0.1.* , PrimitiveArray == 0.8.0.* , PrimitiveArray-Pretty == 0.0.0.* , ShortestPathProblems == 0.0.0.*@@ -79,6 +92,8 @@ BioInf.MutationOrder.EdgeProb BioInf.MutationOrder.MinDist BioInf.MutationOrder.RNA+ BioInf.MutationOrder.SequenceDB+ BioInf.MutationOrder.BackMutations default-extensions: BangPatterns , CPP , DeriveDataTypeable@@ -87,6 +102,7 @@ , GADTs , LambdaCase , MultiParamTypeClasses+ , MultiWayIf , OverloadedStrings , QuasiQuotes , RecordWildCards@@ -95,6 +111,7 @@ , TupleSections , TypeFamilies , TypeOperators+ , UnicodeSyntax default-language: Haskell2010 ghc-options:@@ -111,6 +128,8 @@ build-depends: base , bytestring , cmdargs >= 0.10+ , directory+ , errors , file-embed >= 0.0.8 , filepath --@@ -121,6 +140,8 @@ , DeriveDataTypeable , RecordWildCards , TemplateHaskell+ , UnicodeSyntax+ default-language: main-is: MutationOrder.hs default-language:@@ -165,4 +186,9 @@ source-repository head type: git location: git://github.com/choener/MutationOrder++source-repository this+ type: git+ location: git://github.com/choener/MutationOrder/tree/0.0.1.0+ tag: 0.0.1.0
README.md view
@@ -25,6 +25,23 @@ # Usage instructions +## sequence generation++First, the sequence data base needs to be created. The following assumptions are being made:+- chimp_118.fa is the origin sequence.+- human_118.fa is the target sequence.+- all known mutations are to be ordered.+- One intermediate (or backmutation) is allowed. This will already lead to an+ expansion of the sequence space from ca. 250K sequences to 83.6M sequences!+ Use your local compute cluster or download our precalculated data.++The following command will prepare the working database and populate the seqs subdirectory.++ mkdir workdb+ mkdir workdb/seqs+ mkdir workdb/rnafold+ ./MutationOrder gensequences -w workdb --ancestral chimp_118.fa -e human_118.fa -g 1 --sequencelimit 100000000 --alphabet=ACGT --seqsperfile=100000+ ## example usage We assume that you have two Fasta files, *chimp_118.fa* and *human_118.fa* but@@ -95,10 +112,13 @@ # Installation +Pre-built binaries for Linux are avaiable under [github+releases](https://github.com/choener/MutationOrder/releases)+ 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.+the bottom of the page for instructions to build from source.+ #### Contact
changelog.md view
@@ -1,3 +1,9 @@+0.0.1.0+-------++- new model that includes back-mutations (up to a small constant number)+- compensatory mutations+ 0.0.0.2 -------
src/MutationOrder.hs view
@@ -3,22 +3,19 @@ module Main where +import Control.Monad.IO.Class (liftIO, MonadIO) import Control.Monad+import Control.Error import Data.FileEmbed import qualified Data.ByteString.Char8 as BS import System.Console.CmdArgs+import System.Directory (createDirectoryIfMissing) import System.Exit (exitSuccess, exitFailure) import System.FilePath import BioInf.MutationOrder-import BioInf.MutationOrder.RNA (createRNAlandscape)--data ScoreType- = Mfe- | Centroid- | PairDistMfe- | PairDistCen- deriving (Show,Data,Typeable)+import BioInf.MutationOrder.RNA+import BioInf.MutationOrder.SequenceDB data Options = MutationOrder@@ -37,10 +34,50 @@ , posscaled :: Maybe (Double,Double) , lkupfile :: Maybe FilePath , showmanual :: Bool+ , everykth :: Int } | GenSequences- { infiles :: [FilePath]+ { ancestralSequence ∷ FilePath+ -- ^ the (presumably) ancestral sequence from which to start the mutation+ -- order prediction.+ , extantSequence ∷ FilePath+ -- ^ the extant or target sequence to which to mutate.+ -- + , globalbackmutations ∷ Int+ -- ^ how many global back mutations to allow. For example, if set to @1@,+ -- then for each position in the sequence all "unobserved" nucleotides are+ -- possible with one switch to and one switch back. It is generally a bad+ -- idea to have more than @1@ here and we only allow exactly 1 backmutation+ -- with a specialized algorithm (that still takes quite a while to run).+ , backmutationcolumns ∷ [Int]+ -- ^ Additional columns for backmutations, may overlap with observed+ -- mutations.+ , sequenceLimit ∷ Int+ -- ^ Complain if the number of sequences is above this limit.+ , workdb ∷ FilePath+ -- ^ will write files to @workdb </> seqs@+ , prefixlength ∷ Int+ -- ^ how many characters as prefix for writing+ , seqsperfile ∷ Int+ -- ^ how many sequences to write into each file+ , alphabet ∷ String+ -- ^ the different characters that are allowed }+ | Backmutation+ { ancestralSequence ∷ FilePath+ -- ^ the (presumably) ancestral sequence from which to start the mutation+ -- order prediction.+ , extantSequence ∷ FilePath+ -- ^ the extant or target sequence to which to mutate.+ , workdb ∷ FilePath+ -- ^+ , position ∷ Int+ -- ^+ , alphabet ∷ String+ , scoretype :: ScoreType+ , positivesquared :: Bool+ , onlypositive :: Bool+ } deriving (Show,Data,Typeable) oMutationOrder = MutationOrder@@ -59,26 +96,80 @@ , posscaled = Nothing &= help "--posscaled=x,y scale all values >= x by using y as exponent" , lkupfile = Nothing &= help "developer option: if an RNAfold file with foldings exists, then use it" , showmanual = False &= help "shows the manual"+ , everykth = 1 &= help "for optimal history sampling, choose every kth item only" } oGenSequences = GenSequences- { infiles = def &= args+ { ancestralSequence = def+ , extantSequence = def+ , globalbackmutations = 0+ , backmutationcolumns = []+ , sequenceLimit = 1000000+ , workdb = def+ , prefixlength = 4+ , seqsperfile = 10000+ , alphabet = "" } +oBackmutation = Backmutation+ { ancestralSequence = def+ , extantSequence = def+ , workdb = def+ , position = -1+ , alphabet = ""+ , 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"+ }+ main :: IO () main = do- o <- cmdArgs $ modes [oMutationOrder &= auto, oGenSequences] &= verbosity+ o <- cmdArgs $ modes [oMutationOrder &= auto, oGenSequences, oBackmutation] &= verbosity case o of- MutationOrder{} -> mainProgram o- GenSequences{} -> genSequences o+ MutationOrder{} → mainProgram o+ GenSequences{} → genSequences o+ Backmutation{} → runBackmutation o +-- | This is a simple wrapper around the RNA landscape creation. Landscape+-- creation generates all sequences between ancestral and extant sequence. It+-- will take into account additional global backmutations and further+-- backmutation columns. Note that this can *very easily* lead to combinatorial+-- explosion.+--+-- Needs extra options on: (i) globally active backmutations. For each position+-- try all four nucleotides. (ii) Additional active columns. For each position,+-- try all four nucleotides.+ genSequences o = do let GenSequences{..} = o+ e ← runExceptT $ do+ when (null ancestralSequence) $ throwE "ancestral sequence file?"+ when (null extantSequence) $ throwE "extant sequence file?"+ when (null workdb) $ throwE "work db directory?"+ when (null alphabet) $ throwE "use --alphabet=ACGT (or ACGU if your fasta files are RNA-based)"+ a ← liftIO $ stupidReader ancestralSequence+ e ← liftIO $ stupidReader extantSequence+ (numSeqs, origs, sqs) ← createRNAlandscape2+ alphabet+ (Left $ GlobalBackmutations globalbackmutations)+ (map BackmutationCol backmutationcolumns)+ (Ancestral a) (Extant e)+ unless (numSeqs <= sequenceLimit) $ throwE $ "combinatiorial explosion (" ++ show numSeqs ++ "): reduce search space or allow for higher --sequencelimit"+ -- write out sequences+ let sdb = workdb </> "sequences"+ liftIO $ createDirectoryIfMissing True sdb+ writeSequenceFiles sdb prefixlength seqsperfile+ $ origs ++ concatMap (\(_,_,xs) → xs) sqs+ case e of+ Left err → print err >> exitFailure+ Right () → return ()+ {- ancestral <- stupidReader $ infiles !! 0 current <- stupidReader $ infiles !! 1 let ls = snd $ createRNAlandscape Nothing False ancestral current forM_ ls $ \(k,sq) -> BS.putStrLn sq return ()+ -} embeddedManual = $(embedFile "README.md") @@ -92,22 +183,25 @@ putStrLn "\n\n\nThis program expects exactly two equal-length fasta files as input" exitFailure 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+ runMutationOrder isL fillweight fillstyle scoretype positivesquared posscaled onlypositive cooptcount cooptprint lkupfile outprefix workdb temperature equalStart infiles everykth++runBackmutation ∷ Options → IO ()+runBackmutation Backmutation{..} = do+ e ← runExceptT $ do+ when (null ancestralSequence) $ throwE "ancestral sequence file?"+ when (null extantSequence) $ throwE "extant sequence file?"+ when (null workdb) $ throwE "work db directory?"+ when (null alphabet) $ throwE "use --alphabet=ACGT (or ACGU if your fasta files are RNA-based)"+ a ← liftIO $ Ancestral <$> stupidReader ancestralSequence+ e ← liftIO $ Extant <$> stupidReader extantSequence+ let scaleFun = case scoretype of+ Centroid → centroidDelta' onlypositive positivesquared+ Mfe → mfeDelta' onlypositive positivesquared+ PairDistMfe → mfebpdist' onlypositive positivesquared+ PairDistCen → centroidbpdist' onlypositive positivesquared+ runBackmutationVariants scaleFun workdb alphabet a e position+ case e of+ Left err → print (err ∷ String) >> exitFailure+ Right () → return ()+ return ()