flowsim 0.2.7 → 0.2.8
raw patch · 9 files changed
+106/−81 lines, 9 filesdep ~basedep ~bytestringdep ~cmdargs
Dependency ranges changed: base, bytestring, cmdargs
Files
- flowsim.cabal +34/−21
- src/CloneSim.hs +32/−19
- src/Config.hs +31/−32
- src/EmpFile.hs +2/−2
- src/FlowSim.hs +6/−4
- src/Generations/Empirical.hs too large to diff
- src/Generations/GS20.hs +0/−1
- src/Generations/GenBase.hs +0/−2
- src/HPLCount.hs +1/−0
flowsim.cabal view
@@ -1,35 +1,48 @@ Name: flowsim-Version: 0.2.7+Version: 0.2.8 License: GPL-+Cabal-Version: >= 1.6 Author: Ketil Malde Maintainer: Ketil Malde <ketil@malde.org> Category: Bioinformatics Synopsis: Simulate 454 pyrosequencing-Description: pyrosim - a simulator for 454 pyrosequencing data- Simulate shotgun sequencing with Roche's 454 sequencing technology.+Description: Provides clonesim, a clone simulator that simulates shotgun + genomic clones, and flowsim, that takes the output from clonesim+ (or any other Fasta-formatted file) and generates simulated + reads from them mimicing Roche's 454 pyrosequencing technology,+ writing output in 454's native SFF format. The flowgram+ generation is based on empirical distributions derived from real+ data (although analytic distributions are available too, if you+ prefer). .- The Darcs repository is at <http://malde.org/~ketil/biohaskell/pyrosim>.+ The Darcs repository is at <http://malde.org/~ketil/biohaskell/flowsim>. -HomePage: http://bhlog.malde.org/-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://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 Build-Type: Simple-Tested-with: GHC==6.8.2+Tested-with: GHC==6.8.2, GHC==6.10, GHC==6.12.1 -- Data-files: README-Executable: flowsim-Main-Is: FlowSim.hs-Other-Modules: Statistics, Config, EmpFile- Generations.GenBase, Generations.GS20, - Generations.Titanium, Generations.Empirical-Hs-Source-Dirs: src+Executable flowsim+ 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+ Ghc-Options: -Wall+ Hs-Source-Dirs: src -Executable: hplc-Main-Is: HPLCount.hs-Ghc-Options: -main-is HPLCount-Hs-Source-Dirs: src+Executable hplc+ Build-Depends: bio >= 0.4.4, base >= 4 && < 5, containers, bytestring+ Main-Is: HPLCount.hs+ Ghc-Options: -main-is HPLCount+ Hs-Source-Dirs: src -Executable: clonesim-Main-Is: CloneSim.hs-Hs-Source-Dirs: src+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+ Extensions: DeriveDataTypeable+ Ghc-Options: -Wall+ Hs-Source-Dirs: src
src/CloneSim.hs view
@@ -24,18 +24,18 @@ , input :: [FilePath] } deriving (Data,Typeable,Show) -modes :: Mode Conf-modes = mode $ Conf { - lengths = "Uniform 400 800" &= text "model for clone lengths" & typ "DIST"- , count = 10 &= text "number of reads to generate" & typ "INT"- , input = def &= args & typFile- } &= prog "clonesim"- & text "simulate sequence cloning"- & helpSuffix ["Available distributions (DIST):" - ," Uniform a b, Normal mu sigma, LogNormal mu sigma"]+opts :: Conf+opts = Conf { + lengths = "Uniform 400 800" &= help "model for clone lengths" &= typ "DIST"+ , count = 10 &= help "number of reads to generate" &= typ "INT"+ , input = def &= args &= typFile+ } &= program "clonesim"+ &= summary version+ &= details ["Available distributions (DIST):", + "Uniform a b, Normal mu sigma, LogNormal mu sigma"] version :: String-version = "clonesim v0.2.7, copyright 2010 Ketil Malde" +version = "clonesim v0.2.8, copyright 2010 Ketil Malde" data Dir = Fwd | Rev deriving (Eq,Enum,Bounded) @@ -45,7 +45,7 @@ main :: IO () main = do- conf <- cmdArgs version [modes]+ conf <- cmdArgs opts -- print conf let inf = case input conf of [] -> hPutStrLn stderr "clonesim: reading from stdin" >> hReadFasta stdin@@ -55,26 +55,39 @@ ss <- map (defragSeq . castToNuc) `fmap` inf when (null ss) $ error "No sequences in input, exiting!" let sa = listArray (0,length ss-1) ss- hWriteFasta stdout =<< evalRandIO (simulate conf sa)+ la = listArray (0,length ss) $ scanl (+) 0 $ map seqlength ss+ hWriteFasta stdout =<< evalRandIO (simulate conf sa la) type SeqArray = Array Int (Sequence Nuc)+type LenArray = Array Int Int64 -- | the real 'main'-simulate :: RandomGen g => Conf -> SeqArray -> Rand g [Sequence Nuc]-simulate conf sa = do+simulate :: RandomGen g => Conf -> SeqArray -> LenArray -> Rand g [Sequence Nuc]+simulate conf sa sl = do let ldist = read $ lengths conf- (_,asz) = bounds sa+ (_,alz) = bounds sl+ top = fromIntegral (sl!alz) forM [1..count conf] $ const - (do i <- floor `fmap` sample (Uniform 0 $ 1+fromIntegral asz) -- inclusive?- let s = sa!i- p <- round `fmap` sample (Uniform 0 $ fromIntegral (seqlength s-1))+ (do n <- round `fmap` sample (Uniform 0 top)+ let i = bsearch sl n+ p = n - sl!i+ s = sa!i l <- round `fmap` sample ldist dir <- getRandom return (mkClone dir p l s)) +-- bsearch - returns largest array entry less than search value+bsearch :: LenArray -> Int64 -> Int+bsearch a v = go (bounds a) + where go (x,y) | v >= a!y = y+ | y == x+1 = x+ | a!m <= v = go (m,y)+ | a!m > v = go (x,m)+ where m = (x+y) `div` 2+ mkClone :: Dir -> Int64 -> Int64 -> Sequence Nuc -> Sequence Nuc mkClone dir pos len (Seq h s mq) = - Seq (fromStr (toStr h++" "++show pos++case dir of Fwd -> "+"; Rev -> "-"))+ Seq (fromStr (toStr h++":"++show pos++(case dir of Fwd -> ":fwd"; Rev -> ":rev")++" clonesim")) sd qual where sd = B.take len
src/Config.hs view
@@ -6,7 +6,7 @@ import Data.Char (toLower) import Data.List (isPrefixOf, intersperse)-import System.Console.CmdArgs+import System.Console.CmdArgs as C import System.Directory (doesFileExist) import Control.Monad (when) import System.IO (hPutStrLn, stderr, stdin)@@ -20,7 +20,7 @@ import Generations.Empirical version :: String-version = "flowsim v0.2.7, copyright 2010 Ketil Malde"+version = "flowsim v0.2.8, copyright 2010 Ketil Malde" -- data GenName = Titanium | FLX | GS20 deriving (Typeable, Data, Show, Eq) type GenName = String@@ -30,8 +30,8 @@ , degradation:: Dist , model :: FilePath , qualmethod :: GenName- , discardfilters :: [String]- , trimfilters :: [String]+-- , discardfilters :: [String]+-- , trimfilters :: [String] , flowkey :: String , hplinput :: FilePath , flowlength :: Int@@ -40,32 +40,33 @@ , output :: FilePath } deriving (Typeable, Data, Show, Read) -modes :: Mode Config-modes = mode $ Conf - { generation = "Titanium" &= text "454 generation to simulate" & typ "GEN" & flag "G"- , degradation = def &= text "model for degradation of the flow model" & typ "DIST"- , model = def &= text "empirical distribution for flow generation" & typFile- , qualmethod = def &= text "method for calculating quality" & typ "STRING"- , discardfilters = def &= text "discarding filters to apply" & typ "DFILT"- , trimfilters = def &= text "trimming filters to apply" & typ "TFILT"- , flowkey = def &= text "sequence key to start each read (TCAG)"- , hplinput = def &= text "input genome for HPL count estimate" & typFile- , flowlength = def &= text "number of flow cycles to run"- , flowcycle = def &= text "sequence nucleotides in each flow cycle (TACG)"- , inputs = def &= args & typFile- , output = def &= text "output file"- } &= prog "flowsim"- & text "simulate 454 pyrosequencing."- & helpSuffix hs+opts :: Config +opts = Conf + { generation = "Titanium" &= help "454 generation to simulate" &= typ "GEN" &= C.name "G"+ , degradation = def &= help "model for degradation of the flow model" &= typ "DIST"+ , model = def &= help "empirical distribution for flow generation" &= typFile+ , qualmethod = def &= help "method for calculating quality" &= typ "STRING"+-- , discardfilters = def &= help "discarding filters to apply" &= typ "DFILT"+-- , trimfilters = def &= help "trimming filters to apply" &= typ "TFILT"+ , flowkey = def &= help "sequence key to start each read (TCAG)"+ , hplinput = def &= help "input genome for HPL count estimate" &= typFile+ , flowlength = def &= help "number of flows for each read"+ , flowcycle = def &= help "sequence nucleotides in each flow cycle (TACG)"+ , inputs = def &= args &= typFile+ , output = def &= help "output file"+ } &= program "flowsim"+ &= summary version+ &= details hs where hs = ["Generations (GEN): "++concat (intersperse ", " $ map fst generations) , "Distributions (DIST): Uniform a b, Normal mu sigma, LogNormal mu sigma"- , "Discarding filters (DFILT): ...."- , "Trimming filters (TFILT): ...."]+-- , "Discarding filters (DFILT): ...."+-- , "Trimming filters (TFILT): ...."+ ] mkconf :: IO (Generation,IO [Sequence Nuc],HPLprob,FilePath) mkconf = do- cf <- cmdArgs version [modes]+ cf <- cmdArgs opts let inp = map castToNuc `fmap` if null (inputs cf) then hReadFasta stdin else concat `fmap` mapM readFasta (inputs cf) gen <- mkgen cf whenLoud $ inform ("Generation: "++generation cf ++ "\n" ++ show gen)@@ -95,7 +96,7 @@ setFromStr dflt f str = if null str then f dflt else read str mkgen :: Config -> IO Generation -mkgen (Conf g deg fm qm df tf fk hpl fl fc _inputs _output) = do+mkgen (Conf g deg fm qm {- df tf -} fk _hpl fl fc _inputs _output) = do let gen = lookupGen g setG = setFromGen gen setD = setFromStr gen @@ -119,10 +120,10 @@ { qcall = my_qcall , degrade = setD degrade deg , models = my_model- , disc_filters = setF df - , trim_filters = setT tf+-- , disc_filters = setF df +-- , trim_filters = setT tf , f_key = setG f_key fk- , f_len = if fl > 0 then fl else f_len gen+ , f_len = if fl > 0 then 4*((fl+3)`div` 4) else f_len gen -- is this sufficient? , f_cycle = setG f_cycle fc } @@ -132,16 +133,14 @@ -- filterlist :: [DiscardFilter] filterlist = [] -- todo: export from SFF_filters --- trimlist :: [Trimfilter]+-- trimlist :: [TrimFilter] trimlist = [] -- Logging inform :: String -> IO () inform = hPutStrLn stderr -whenLoud, whenNotQuiet :: IO () -> IO ()-whenLoud action = flip when action =<< isLoud--- note that isQuiet = return False in CmdArgs!+whenNotQuiet :: IO () -> IO () whenNotQuiet action = do n <- isNormal l <- isLoud when (n || l) action
src/EmpFile.hs view
@@ -1,7 +1,6 @@ module EmpFile where import Statistics-import Control.Monad (when) readPdf :: FilePath -> IO [Distribution] readPdf f = do @@ -11,10 +10,11 @@ parsePdf :: String -> [String] -> [Distribution] parsePdf f cs = let (_:fs:vss) = transpose (map words cs)- in if (or $ zipWith (/=) [0..] (map read fs))+ in if (or $ zipWith (/=) [(0::Int)..] (map read fs)) then error ("Flow index from "++f++" isn't [0..]!") else map (fromPdf 0 0.01 . map (P . read)) vss +transpose :: [[a]] -> [[a]] transpose xs = let r0 = map head xs rest = map tail xs in if any null rest then [r0] else r0 : transpose rest
src/FlowSim.hs view
@@ -72,9 +72,10 @@ makeModels :: RandomGen g => Generation -> Rand g [Model] makeModels gen = case models gen of [m] -> iterativeWorsen (f_len gen) (degrade gen) m- im@(_:_) -> do- let ns = f_len gen `div` length im- return $ concat $ map (replicate ns) im+ im@(i1:is) -> do+ let (ns,rs) = f_len gen `divMod` length im+ -- error (show (ns,rs,length (replicate (ns+rs) i1 ++ concatMap (replicate ns) is)))+ return (replicate (ns+rs) i1 ++ concatMap (replicate ns) is) iterativeWorsen :: RandomGen g => Int -> Distribution -> Model -> Rand g [Model] iterativeWorsen 0 _ _ = return []@@ -110,6 +111,7 @@ permuteAndCall gen hplc fs = do -- when (trace (show fs) False) (return ()) let n = fromIntegral $ f_len gen+ -- tack on the correct number of zero flows after end of clone flow_cont = drop (length fs `mod` length (f_cycle gen)) $ cycle (f_cycle gen) xs = fs ++ zip flow_cont (replicate (n-length fs) 0) ms <- makeModels gen@@ -173,7 +175,7 @@ -- | Consistency check on generated ReadBlocks. verifyRB :: Int16 -> ReadBlock -> ReadBlock verifyRB fl rb - | name_length rh == 0 = err "name_length is zero"+ -- | 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"
src/Generations/Empirical.hs view
file too large to diff
src/Generations/GS20.hs view
@@ -5,7 +5,6 @@ import Data.Array import qualified Data.ByteString.Char8 as BC import qualified Data.ByteString as B-import Control.Monad (when) type QualTable = Array Flow B.ByteString
src/Generations/GenBase.hs view
@@ -11,8 +11,6 @@ import Bio.Sequence.SFF import Bio.Sequence.SFF_filters -import Debug.Trace (trace)- import Statistics import HPLCount hiding (main)
src/HPLCount.hs view
@@ -17,6 +17,7 @@ mkHPLprobs ft = let totals = [(c, fromIntegral $ sum [M.findWithDefault 0 (c,x) ft | x <- [0..50]]) | c <- "ACGT" ] in \ch n -> fromIntegral (M.findWithDefault 0 (ch,n) ft) / fromJust (lookup ch totals) +main :: IO () main = do [f] <- getArgs display "tacg" . freqtable . concatMap (hpls "tacg" . seqdata) =<< readFasta f