diff --git a/flowsim.cabal b/flowsim.cabal --- a/flowsim.cabal +++ b/flowsim.cabal @@ -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 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 . + The Darcs repository is at . -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 diff --git a/src/CloneSim.hs b/src/CloneSim.hs --- a/src/CloneSim.hs +++ b/src/CloneSim.hs @@ -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 diff --git a/src/Config.hs b/src/Config.hs --- a/src/Config.hs +++ b/src/Config.hs @@ -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 diff --git a/src/EmpFile.hs b/src/EmpFile.hs --- a/src/EmpFile.hs +++ b/src/EmpFile.hs @@ -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 diff --git a/src/FlowSim.hs b/src/FlowSim.hs --- a/src/FlowSim.hs +++ b/src/FlowSim.hs @@ -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" diff --git a/src/Generations/Empirical.hs b/src/Generations/Empirical.hs # file too large to diff: src/Generations/Empirical.hs diff --git a/src/Generations/GS20.hs b/src/Generations/GS20.hs --- a/src/Generations/GS20.hs +++ b/src/Generations/GS20.hs @@ -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 diff --git a/src/Generations/GenBase.hs b/src/Generations/GenBase.hs --- a/src/Generations/GenBase.hs +++ b/src/Generations/GenBase.hs @@ -11,8 +11,6 @@ import Bio.Sequence.SFF import Bio.Sequence.SFF_filters -import Debug.Trace (trace) - import Statistics import HPLCount hiding (main) diff --git a/src/HPLCount.hs b/src/HPLCount.hs --- a/src/HPLCount.hs +++ b/src/HPLCount.hs @@ -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