flowsim 0.2.8 → 0.3
raw patch · 15 files changed
+243/−49 lines, 15 filesdep ~bionew-component:exe:duplicatornew-component:exe:gelfilternew-component:exe:kitsimnew-component:exe:mutator
Dependency ranges changed: bio
Files
- flowsim.cabal +34/−6
- src/CloneSim.hs +2/−4
- src/Config.hs +2/−4
- src/Duplicator.hs +24/−0
- src/FlowSim.hs +13/−14
- src/GelFilter.hs +14/−0
- src/Generations/Empirical.hs too large to diff
- src/Generations/GS20.hs +3/−2
- src/Generations/GenBase.hs +10/−9
- src/Generations/Titanium.hs +18/−2
- src/HPLCount.hs +1/−0
- src/KitSim.hs +43/−0
- src/Mutator.hs +68/−0
- src/Statistics.hs +7/−8
- src/Version.hs +4/−0
flowsim.cabal view
@@ -1,5 +1,5 @@ Name: flowsim-Version: 0.2.8+Version: 0.3 License: GPL Cabal-Version: >= 1.6 Author: Ketil Malde@@ -18,8 +18,7 @@ . The Darcs repository is at <http://malde.org/~ketil/biohaskell/flowsim>. -HomePage: http://blog.malde.org/index.php/flowsim/--- Build-Depends: bio >= 0.4.4, base >= 4 && <= 5, array >= 0.1, bytestring >= 0.9.1, random, MonadRandom, cmdargs <= 0.1.1, containers, directory+HomePage: http://biohaskell.org/Applications/FlowSim Build-Type: Simple Tested-with: GHC==6.8.2, GHC==6.10, GHC==6.12.1 @@ -28,8 +27,8 @@ Main-Is: FlowSim.hs Other-Modules: Statistics, Config, EmpFile Generations.GenBase, Generations.GS20, - Generations.Titanium, Generations.Empirical- Build-Depends: bio >= 0.4.4, base >= 4 && < 5, array >= 0.1, random, MonadRandom, cmdargs >= 0.5, containers, bytestring, directory+ Generations.Titanium, Generations.Empirical, Version+ Build-Depends: bio >= 0.4.9, base >= 4 && < 5, array >= 0.1, random, MonadRandom, cmdargs >= 0.5, containers, bytestring, directory Ghc-Options: -Wall Hs-Source-Dirs: src @@ -42,7 +41,36 @@ Executable clonesim Main-Is: CloneSim.hs Build-Depends: bio >= 0.4.4, base >= 4 && < 5, array >= 0.1, cmdargs >= 0.5, bytestring, MonadRandom- Other-Modules: Statistics+ Other-Modules: Statistics, Version Extensions: DeriveDataTypeable Ghc-Options: -Wall+ Hs-Source-Dirs: src++Executable kitsim+ Main-Is: KitSim.hs+ Build-Depends: bio >= 0.4.9, base >= 4 && < 5, cmdargs >= 0.5+ Other-Modules: Version+ Extensions: DeriveDataTypeable+ Ghc-Options: -Wall -main-is KitSim+ Hs-Source-Dirs: src++Executable mutator+ Main-Is: Mutator.hs+ Build-Depends: bio >= 0.4.9, base >= 4 && < 5, cmdargs >= 0.5+ Other-Modules: Version, Statistics+ Extensions: DeriveDataTypeable+ Ghc-Options: -Wall -main-is Mutator+ Hs-Source-Dirs: src++Executable duplicator+ Main-Is: Duplicator.hs+ Build-Depends: bio >= 0.4.9, base >= 4 && < 5+ Other-Modules: Statistics+ Ghc-Options: -Wall -main-is Duplicator+ Hs-Source-Dirs: src++Executable gelfilter+ Main-Is: GelFilter.hs+ Build-Depends: bio >= 0.4.9, base >= 4 && < 5+ Ghc-Options: -Wall -main-is GelFilter Hs-Source-Dirs: src
src/CloneSim.hs view
@@ -9,6 +9,7 @@ import Bio.Sequence hiding ((!)) import Statistics+import Version import System.IO (stdin,stdout,stderr,hPutStrLn) import Control.Monad (forM,when)@@ -30,12 +31,9 @@ , count = 10 &= help "number of reads to generate" &= typ "INT" , input = def &= args &= typFile } &= program "clonesim"- &= summary version+ &= summary ("clonesim "++version) &= details ["Available distributions (DIST):", "Uniform a b, Normal mu sigma, LogNormal mu sigma"]--version :: String-version = "clonesim v0.2.8, copyright 2010 Ketil Malde" data Dir = Fwd | Rev deriving (Eq,Enum,Bounded)
src/Config.hs view
@@ -11,6 +11,7 @@ import Control.Monad (when) import System.IO (hPutStrLn, stderr, stdin) import EmpFile+import Version import Bio.Sequence @@ -19,9 +20,6 @@ import Generations.Titanium import Generations.Empirical -version :: String-version = "flowsim v0.2.8, copyright 2010 Ketil Malde"- -- data GenName = Titanium | FLX | GS20 deriving (Typeable, Data, Show, Eq) type GenName = String type Dist = String@@ -55,7 +53,7 @@ , inputs = def &= args &= typFile , output = def &= help "output file" } &= program "flowsim"- &= summary version+ &= summary ("flowsim "++version) &= details hs where hs = ["Generations (GEN): "++concat (intersperse ", " $ map fst generations)
+ src/Duplicator.hs view
@@ -0,0 +1,24 @@+{- Extremely simple program to introduce duplicates -}++module Duplicator where++import System.IO+import System.Environment (getArgs)++import Bio.Sequence+import Statistics++main :: IO ()+main = do+ args <- getArgs+ case args of [p] -> hReadFasta stdin >>= evalRandIO . dup (read p) . map castToNuc >>= hWriteFasta stdout+ _ -> error "Usage: duplicator pr\n where pr is the probability of recursively duplicating a sequence"++dup :: RandomGen g => Double -> [Sequence Nuc] -> Rand g [Sequence Nuc]+dup p (x:xs) = do+ r <- getRandomR (0,1)+ if r < p + then dup p (x:x:xs)+ else do ys <- dup p xs+ return (x:ys)+dup _ [] = return []
src/FlowSim.hs view
@@ -12,7 +12,6 @@ - sort and generate -}-{-# options -XParallelListComp #-} module Main where @@ -47,7 +46,7 @@ (gen,is,hplc,o) <- mkconf ss <- is- case ss of [] -> error ("Input appears to be empty?")+ case ss of [] -> error "Input appears to be empty?" _ -> return () sff <- evalRandIO $ sim454 gen hplc ss n <- writeSFF' o sff@@ -61,9 +60,9 @@ sim454 :: RandomGen g => Generation -> HPLprob -> [Sequence Nuc] -> Rand g SFF sim454 gen hplc ss = do let ch = makeCommonHeader gen- tf r = foldr ($) r $ trim_filters gen - df r = and $ zipWith ($) (disc_filters gen) (repeat r)- rbs <- map tf `fmap` filter df `fmap` mapM (\s -> makeReadBlock gen hplc s ch) ss+ tf r = foldr ($) r (trim_primer (adapter gen) : trim_filters gen)+ df r = and $ zipWith ($) (discard_key (f_key gen) : disc_filters gen) (repeat r)+ rbs <- map tf `fmap` filter df `fmap` mapM (makeReadBlock gen hplc ch) ss return (SFF ch rbs) -- Generate a sequence of models@@ -96,7 +95,7 @@ makeFlows :: [Char] -> SeqData -> [(Char,Int)] makeFlows c s | B.null s = [] | otherwise = let (c1,s1) = makeCycle [] (take 4 c) s - s2 = case B.uncons s1 of Just (x,_) -> if not (toUpper x `elem` "ACGT") then B.tail s1 else s1; _ -> s1+ s2 = case B.uncons s1 of Just (x,_) -> if notElem (toUpper x) "ACGT" then B.tail s1 else s1; _ -> s1 in c1 ++ makeFlows c s2 makeCycle :: [(Char,Int)] -> [Char] -> SeqData -> ([(Char,Int)],SeqData)@@ -149,13 +148,13 @@ -- | Generate a ReadBlock -- direction and position chosen at random, and encoded in the read name-makeReadBlock :: RandomGen g => Generation -> HPLprob -> Sequence Nuc -> CommonHeader -> Rand g ReadBlock-makeReadBlock g hplc sq ch = do+makeReadBlock :: RandomGen g => Generation -> HPLprob -> CommonHeader -> Sequence Nuc -> Rand g ReadBlock+makeReadBlock g hplc ch sq = do let sdata = seqdata sq rn = BS.concat $ BL.toChunks $ seqlabel sq- fs = makeFlows (BS.unpack $ flow ch) $ BL.concat [BL.fromChunks [key ch], {- B.take (floor cl) $ -} sdata]+ fs = makeFlows (BS.unpack $ flow ch) sdata (pfs,cs,qs,is) <- convertCalls `fmap` permuteAndCall g hplc fs- return $ verifyRB (flow_length ch) $ ReadBlock {+ return $ verifyRB (flow_length ch) ReadBlock { read_header = ReadHeader { name_length = fromIntegral $ BS.length rn -- :: Int16 , num_bases = fromIntegral $ length cs -- :: Int32@@ -178,11 +177,11 @@ -- | name_length rh == 0 = err "name_length is zero" | num_bases rh == 0 = err "num_bases is zero" -- | clip_qual_left rh > clip_qual_right rh = err "clipping gives negative sequence"- | BS2.length (read_name rh) /= (fromIntegral $ name_length rh) = err "read_name has incorrect length"+ | BS2.length (read_name rh) /= fromIntegral (name_length rh) = err "read_name has incorrect length" | BS.length (flow_data rb) `div` 2 /= fromIntegral fl = err ("Number of flows ("++show (BS.length (flow_data rb) `div` 2)++")do not match flow_length of "++show fl++" in CommonHeader")- | BS.length (flow_data rb) `div` 2 < (fromIntegral $ sum $ BS2.unpack $ flow_index rb) = err "flow_index longer than flows"- | B.length (bases rb) /= (fromIntegral $ BS.length $ flow_index rb) = err ("bases ("++show (B.length (bases rb))++") and flow_index ("++show (BS.length $ flow_index rb)++") have different lengths")- | B.length (quality rb) /= (fromIntegral $ BS.length $ flow_index rb) = err ("quality ("++show (B.length (quality rb))++") and flow_index ("++show (BS.length $ flow_index rb)++") have different lengths")+ | BS.length (flow_data rb) `div` 2 < fromIntegral (sum $ BS2.unpack $ flow_index rb) = err "flow_index longer than flows"+ | B.length (bases rb) /= fromIntegral (BS.length $ flow_index rb) = err ("bases ("++show (B.length (bases rb))++") and flow_index ("++show (BS.length $ flow_index rb)++") have different lengths")+ | B.length (quality rb) /= fromIntegral (BS.length $ flow_index rb) = err ("quality ("++show (B.length (quality rb))++") and flow_index ("++show (BS.length $ flow_index rb)++") have different lengths") | otherwise = rb where rh = read_header rb err str = error (str ++ "\n" ++ show rb)
+ src/GelFilter.hs view
@@ -0,0 +1,14 @@+{- Extremely simple tool to filter reads by length -}++module GelFilter where++import System.Environment (getArgs)+import System.IO++import Bio.Sequence++main :: IO ()+main = do+ [min,max] <- map read `fmap` getArgs+ hWriteFasta stdout =<< filter (\x -> seqlength x>=min && seqlength x<= max) `fmap` hReadFasta stdin+
src/Generations/Empirical.hs view
file too large to diff
src/Generations/GS20.hs view
@@ -29,12 +29,13 @@ in if h == 0 then LogNormal (negate 2.5) 0.2 else Normal h' (0.04+0.06*h')] , qcall = qual_gs20_tab , degrade = Normal 0.005 0.002- , disc_filters = [filter_key,filter_empty]- , trim_filters = []+ , disc_filters = [discard_empty]+ , trim_filters = [trim_sigint, trim_qual20 10] , f_key = "TCAG" , f_len = 168 -- :: Int16 , f_cycle = "TACG"+ , adapter = "ctgactgagacacgcaacaggggataggacaaggcacacagggggatagg" -- taken from DZX0XNV01.sff }
src/Generations/GenBase.hs view
@@ -30,6 +30,7 @@ , f_len :: Int , f_cycle :: String -- readblock :: [ReadBlock] ?+ , adapter :: String } instance Show Generation where@@ -43,21 +44,21 @@ -- | Calculate exact qualities from model using Bayes' theorem -- This should emulate Marguiles et al.'s method. qual_exact_decreasing :: HPLprob -> Model -> (Char,Flow) -> [Qual]-qual_exact_decreasing ft m (c,flow) = - let f = fromIntegral flow / 100- prob hpl | hpl > 20 = error "can't handle hpls > 20" - | otherwise = ft c hpl * pdf (m hpl) f / sum [pdf (m x) f * ft c x | x <- [0..20]]+qual_exact_decreasing ft m (c,fl) =+ let f = fromIntegral fl / 100+ h = fromIntegral $ (fl+50) `div` 100+ prob hpl = ft c hpl * pdf (m hpl) f / sum [pdf (m x) f * ft c x | x <- [0..20]] probs = drop 1 $ scanl (+) 0 $ map prob [0..20]- quals = map (round . max 0 . min 60 . ((-10)*) . logBase 10) $ probs+ quals = replicate (h-20) 60 ++ (map (round . max 0 . min 60 . ((-10)*) . logBase 10) $ probs) in -- trace ("#"++show (c,flow) ++"\n"++show probs ++"\n"++show quals) $- take (fromIntegral $ (flow+50) `div` 100) $ quals+ take (fromIntegral $ (fl+50) `div` 100) $ quals -- | Using Bayes, but only calculating the probablity of the call lenght -- being correct. Should be similar to Titanium (but is it?) qual_exact_fixed :: HPLprob -> Model -> (Char,Flow) -> [Qual]-qual_exact_fixed ft m (c,flow) = - let f = fromIntegral flow / 100- hpl = fromIntegral $ (flow+50) `div` 100+qual_exact_fixed ft m (c,fl) = + let f = fromIntegral fl / 100+ hpl = fromIntegral $ (fl+50) `div` 100 prob = ft c hpl * pdf (m hpl) f / sum [pdf (m x) f * ft c x | x <- [0..20]] in replicate hpl (round $ max 0 . min 60 $ (-10) * logBase 10 (1-prob))
src/Generations/Titanium.hs view
@@ -1,5 +1,7 @@ module Generations.Titanium where+ import Generations.GenBase+import qualified Data.ByteString.Char8 as BC titanium :: Generation titanium = Gen@@ -8,9 +10,23 @@ , models = [\h -> let h' = fromIntegral h in if h == 0 then LogNormal (negate 2.5) 0.2 else Normal h' (0.04+0.05*h')] , degrade = Normal 0.001 0.001- , disc_filters = [filter_key,filter_empty]- , trim_filters = []+ , disc_filters = [discard_empty,discard_dots 0.05,discard_mixed,discard_length 186] -- discard_key+ , trim_filters = [trim_sigint, trim_qual20 10] -- find_primer , f_key = "TCAG" , f_len = 800 -- :: Int16 , f_cycle = "TACG"+ , adapter = "ctgagactgccaaggcacacagggggatagg" }++{- seen:+ ggcgggcgatgtctcgtctgagcgggctggc aaggc <- in transcriptome seq.+ ctgagactgcc aaggcacacagggggatagg <- sea bass and cod?+ agtcgtggaggc aaggcacacaggg--atagg <- "rapid" kit++From: http://seqanswers.com/forums/showthread.php?t=2698+ TCGTATAACTTCGTATAATGTATGCTATACGAAGTTATTACG + GTTGGAACCGAAAGGGTTTGAATTCAAACCCTTTCGGTTCCAAC++Todo: check MIRA and sffToWhatever from Celera.++-}
src/HPLCount.hs view
@@ -43,3 +43,4 @@ | not (elem (toUpper $ B.head sd) xs') = go (f:fs) (B.tail sd) | otherwise = let (this,rest) = B.span ((==f).toUpper) sd in (f,fromIntegral (B.length this)) : go fs rest+ go [] _ = error "should never happen"
+ src/KitSim.hs view
@@ -0,0 +1,43 @@+{-| Simulate reagent 'kits', i.e. tack on adapters etc -}+{-# LANGUAGE DeriveDataTypeable #-}++module KitSim where++import Version++import System.Console.CmdArgs+import System.IO+import qualified Data.ByteString.Lazy.Char8 as B++import Bio.Sequence++main :: IO ()+main = do+ cf <- cmdArgs conf+ let reader = fmap (map castToNuc) $+ case input cf of "-" -> hReadFasta stdin+ f -> readFasta f+ writer = case output cf of "-" -> hWriteFasta stdout+ f -> writeFasta f+ apply cf `fmap` reader >>= writer++apply :: Conf -> [Sequence Nuc] -> [Sequence Nuc]+apply c = map apply1+ where apply1 (Seq s d _) = Seq s (B.concat [k,d,a]) Nothing+ k = B.pack (key c)+ a = B.pack (adapter c)++conf :: Conf+conf = Conf + { key = "TCAG" &= help "Sequence for initial key" &= typ "String"+ , adapter = "ctgagactgccaaggcacacagggggatagg" &= help "Sequence for B-adapter" &= typ "String"+ , input = "-" &= args &= typFile+ , output = "-" &= help "Output file" &= typFile+ } &= program "kitsim"+ &= summary ("kitsim"++version)+ &= details ["Simulates the sequencing kit by tacking on the initial key"+ ,"(really the end of the A-adapter) and B-adapter"]++data Conf = Conf { key, adapter :: String, input, output :: FilePath }+ deriving (Data,Typeable)+
+ src/Mutator.hs view
@@ -0,0 +1,68 @@+{-| The mutator randomly introduces substitutions and indels into + Fasta sequences+-}+{-# LANGUAGE DeriveDataTypeable #-}++module Mutator where++import System.IO+import System.Console.CmdArgs+import Bio.Sequence++import Statistics+import Version++data Conf = Conf { input, output :: FilePath+ , subst, indel :: Double+ -- , gapext :: Double+ } deriving (Typeable, Data)++conf :: Conf+conf = Conf + { input = "-" &= args &= typFile+ , output = "-" &= help "Output file" &= typFile+ , subst = 0.01 &= help "Substitution rate" &= typ "Float"+ , indel = 0.01 &= help "Indel rate" &= typ "Float"+ -- , gapext = 0.05 &= help "Gap extension rate" &= typ "Float" -- todo: affine gaps+ } &= program "mutator"+ &= summary ("mutator "++version)+ &= details ["Mutate sequences in Fasta format by introducing"+ ,"random substitutions and insertions/deletions"]++main :: IO ()+main = do+ c <- cmdArgs conf+ let inp = case input c of "-" -> hReadFasta stdin+ x -> readFasta x+ outp = case output c of "-" -> hWriteFasta stdout+ x -> writeFasta x+ inp >>= doMutate c . map castToNuc >>= outp+ +doMutate :: Conf -> [Sequence Nuc] -> IO [Sequence Nuc]+doMutate cf = evalRandIO . mutate (subst cf) (indel cf)++mutate :: RandomGen g => + Double -> Double -> [Sequence Nuc] -> Rand g [Sequence Nuc]+mutate sub ind = mapM mut1 + where mut1 (Seq h d _) = do + d2 <- go (toStr d)+ return $ Seq h (fromStr d2) Nothing+ nuc :: Int -> Char + nuc x = case x of {0 -> 'a'; 1 -> 'c'; 2 -> 'g'; 3 -> 't'}+ go "" = return ""+ go (x:xs) = do+ z <- sample (Uniform 0 1)+ if z<sub + then do -- substitute+ y <- getRandomR (0,3)+ ys <- go xs+ return (nuc y:ys)+ else if z < sub+(ind/2) then do -- insert+ y <- getRandomR (0,3)+ ys <- go (x:xs)+ return (nuc y:ys)+ else if z < sub+ind then do -- delete+ go xs+ else do+ ys <- go xs+ return (x:ys)
src/Statistics.hs view
@@ -35,7 +35,7 @@ -- to uniformly spaced points starting at 'start' with 'step' points per unit. -- Automatically center on pvalue = 50. fromPdf :: Double -> Double -> [Prob] -> Distribution-fromPdf start h ps = let ps' = acc 0 $ map ((/scale) . unprob) $ ps+fromPdf start h ps = let ps' = acc 0 $ map ((/scale) . unprob) ps mu = start+(fromIntegral . length . takeWhile (<0.5)) ps'*h-h/2 a = floor ((start-mu)/h) b = a+length ps'@@ -96,13 +96,12 @@ | otherwise = 0 pdf (Empirical mu h cd) x = let (a,b) = bounds cd x' = (x-mu)/h- in if x' < fromIntegral a then 0- else if x' >= fromIntegral b then 0- else let x1 = floor x'- x2 = x1+1- y1 = cd!x1- y2 = cd!x2- in (y2-y1)/h+ in if x' < fromIntegral a || x' >= fromIntegral b then 0+ else let x1 = floor x'+ x2 = x1+1+ y1 = cd!x1+ y2 = cd!x2+ in (y2-y1)/h -- ------------------------------ -- Specifics -- ------------------------------
+ src/Version.hs view
@@ -0,0 +1,4 @@+module Version where++version :: String+version = " v0.3, copyright 2010 Ketil Malde"