diff --git a/flowsim.cabal b/flowsim.cabal --- a/flowsim.cabal +++ b/flowsim.cabal @@ -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 . -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 diff --git a/src/CloneSim.hs b/src/CloneSim.hs --- a/src/CloneSim.hs +++ b/src/CloneSim.hs @@ -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) diff --git a/src/Config.hs b/src/Config.hs --- a/src/Config.hs +++ b/src/Config.hs @@ -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) diff --git a/src/Duplicator.hs b/src/Duplicator.hs new file mode 100644 --- /dev/null +++ b/src/Duplicator.hs @@ -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 [] diff --git a/src/FlowSim.hs b/src/FlowSim.hs --- a/src/FlowSim.hs +++ b/src/FlowSim.hs @@ -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) diff --git a/src/GelFilter.hs b/src/GelFilter.hs new file mode 100644 --- /dev/null +++ b/src/GelFilter.hs @@ -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 + 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 @@ -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 } diff --git a/src/Generations/GenBase.hs b/src/Generations/GenBase.hs --- a/src/Generations/GenBase.hs +++ b/src/Generations/GenBase.hs @@ -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)) diff --git a/src/Generations/Titanium.hs b/src/Generations/Titanium.hs --- a/src/Generations/Titanium.hs +++ b/src/Generations/Titanium.hs @@ -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. + +-} diff --git a/src/HPLCount.hs b/src/HPLCount.hs --- a/src/HPLCount.hs +++ b/src/HPLCount.hs @@ -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" diff --git a/src/KitSim.hs b/src/KitSim.hs new file mode 100644 --- /dev/null +++ b/src/KitSim.hs @@ -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) + diff --git a/src/Mutator.hs b/src/Mutator.hs new file mode 100644 --- /dev/null +++ b/src/Mutator.hs @@ -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 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 -- ------------------------------ diff --git a/src/Version.hs b/src/Version.hs new file mode 100644 --- /dev/null +++ b/src/Version.hs @@ -0,0 +1,4 @@ +module Version where + +version :: String +version = " v0.3, copyright 2010 Ketil Malde"