diff --git a/BioInf/MutationOrder.hs b/BioInf/MutationOrder.hs --- a/BioInf/MutationOrder.hs +++ b/BioInf/MutationOrder.hs @@ -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 () +-} diff --git a/BioInf/MutationOrder/BackMutations.hs b/BioInf/MutationOrder/BackMutations.hs new file mode 100644 --- /dev/null +++ b/BioInf/MutationOrder/BackMutations.hs @@ -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 + + + diff --git a/BioInf/MutationOrder/MinDist.hs b/BioInf/MutationOrder/MinDist.hs --- a/BioInf/MutationOrder/MinDist.hs +++ b/BioInf/MutationOrder/MinDist.hs @@ -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 = diff --git a/BioInf/MutationOrder/RNA.hs b/BioInf/MutationOrder/RNA.hs --- a/BioInf/MutationOrder/RNA.hs +++ b/BioInf/MutationOrder/RNA.hs @@ -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) -- | -- diff --git a/BioInf/MutationOrder/SequenceDB.hs b/BioInf/MutationOrder/SequenceDB.hs new file mode 100644 --- /dev/null +++ b/BioInf/MutationOrder/SequenceDB.hs @@ -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 #-} + diff --git a/MutationOrder.cabal b/MutationOrder.cabal --- a/MutationOrder.cabal +++ b/MutationOrder.cabal @@ -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 diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -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 diff --git a/changelog.md b/changelog.md --- a/changelog.md +++ b/changelog.md @@ -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 ------- diff --git a/src/MutationOrder.hs b/src/MutationOrder.hs --- a/src/MutationOrder.hs +++ b/src/MutationOrder.hs @@ -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 ()