diff --git a/flowsim.cabal b/flowsim.cabal --- a/flowsim.cabal +++ b/flowsim.cabal @@ -1,6 +1,6 @@ Name: flowsim -Version: 0.3 -License: GPL +Version: 0.3.3 +License: GPL-2 Cabal-Version: >= 1.6 Author: Ketil Malde Maintainer: Ketil Malde @@ -20,7 +20,7 @@ HomePage: http://biohaskell.org/Applications/FlowSim Build-Type: Simple -Tested-with: GHC==6.8.2, GHC==6.10, GHC==6.12.1 +Tested-with: GHC==6.8.2, GHC==6.10, GHC==6.12.1, GHC==7.4.1 -- Data-files: README Executable flowsim @@ -28,19 +28,19 @@ Other-Modules: Statistics, Config, EmpFile Generations.GenBase, Generations.GS20, 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 + Build-Depends: base >= 4 && < 5, array >= 0.1, random, MonadRandom, cmdargs >= 0.5, containers, bytestring, directory, mtl >= 2, biofasta, biocore >= 0.3, biosff >= 0.3 Ghc-Options: -Wall Hs-Source-Dirs: src Executable hplc - Build-Depends: bio >= 0.4.4, base >= 4 && < 5, containers, bytestring + Build-Depends: base >= 4 && < 5, containers, bytestring Main-Is: HPLCount.hs Ghc-Options: -main-is HPLCount 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 + Build-Depends: base >= 4 && < 5, array >= 0.1, cmdargs >= 0.5, bytestring, MonadRandom Other-Modules: Statistics, Version Extensions: DeriveDataTypeable Ghc-Options: -Wall @@ -48,7 +48,7 @@ Executable kitsim Main-Is: KitSim.hs - Build-Depends: bio >= 0.4.9, base >= 4 && < 5, cmdargs >= 0.5 + Build-Depends: base >= 4 && < 5, cmdargs >= 0.5 Other-Modules: Version Extensions: DeriveDataTypeable Ghc-Options: -Wall -main-is KitSim @@ -56,7 +56,7 @@ Executable mutator Main-Is: Mutator.hs - Build-Depends: bio >= 0.4.9, base >= 4 && < 5, cmdargs >= 0.5 + Build-Depends: base >= 4 && < 5, cmdargs >= 0.5 Other-Modules: Version, Statistics Extensions: DeriveDataTypeable Ghc-Options: -Wall -main-is Mutator @@ -64,13 +64,19 @@ Executable duplicator Main-Is: Duplicator.hs - Build-Depends: bio >= 0.4.9, base >= 4 && < 5 + Build-Depends: 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 + Build-Depends: base >= 4 && < 5 Ghc-Options: -Wall -main-is GelFilter + Hs-Source-Dirs: src + +Executable filtersff + Main-Is: FilterSFF.hs + Build-Depends: base >= 4 && < 5 + Ghc-Options: -Wall -main-is FilterSFF Hs-Source-Dirs: src diff --git a/src/CloneSim.hs b/src/CloneSim.hs --- a/src/CloneSim.hs +++ b/src/CloneSim.hs @@ -7,16 +7,18 @@ module Main where -import Bio.Sequence hiding ((!)) +import Bio.Core.Sequence +import Bio.Sequence.Fasta + import Statistics import Version import System.IO (stdin,stdout,stderr,hPutStrLn) import Control.Monad (forM,when) import qualified Data.ByteString.Lazy as B -import Data.Int (Int64) import System.Console.CmdArgs import Data.Array +import Data.Char type Dist = String data Conf = Conf @@ -50,32 +52,32 @@ [x] -> readFasta x _ -> error "Please specify only a single input file" seq inf $ return () -- force the above warning - ss <- map (defragSeq . castToNuc) `fmap` inf + ss <- map defragSeq `fmap` inf when (null ss) $ error "No sequences in input, exiting!" let sa = listArray (0,length ss-1) ss 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 +type SeqArray = Array Int Sequence +type LenArray = Array Int Offset + -- | the real 'main' -simulate :: RandomGen g => Conf -> SeqArray -> LenArray -> Rand g [Sequence Nuc] +simulate :: RandomGen g => Conf -> SeqArray -> LenArray -> Rand g [Sequence] simulate conf sa sl = do let ldist = read $ lengths conf (_,alz) = bounds sl top = fromIntegral (sl!alz) forM [1..count conf] $ const - (do n <- round `fmap` sample (Uniform 0 top) + (do n <- floor `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)) + return $ mkClone dir p (Offset l) s) -- bsearch - returns largest array entry less than search value -bsearch :: LenArray -> Int64 -> Int +bsearch :: LenArray -> Offset -> Int bsearch a v = go (bounds a) where go (x,y) | v >= a!y = y | y == x+1 = x @@ -83,13 +85,31 @@ | 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 -> ":fwd"; Rev -> ":rev")++" clonesim")) - sd - qual - where sd = B.take len +mkClone :: Dir -> Offset -> Offset -> Sequence -> Sequence +mkClone dir (Offset pos) (Offset len) (Seq h (SeqData s) mq) = + Seq label (SeqData sd) (qual mq) + where label = fromString (toString h++":"++show pos++ + (case dir of Fwd -> ":fwd"; Rev -> ":rev")++" clonesim") + sd = B.take len . (case dir of Fwd -> snd; Rev -> revcompl' . fst) . B.splitAt pos $ s - qual = B.take len `fmap` case dir of Fwd -> fst; Rev -> snd - `fmap` B.splitAt pos `fmap` mq + qual Nothing = Nothing + qual (Just (QualData qd)) = Just . QualData . B.take len + . case dir of Fwd -> fst; Rev -> snd + $ B.splitAt pos qd + +revcompl' :: B.ByteString -> B.ByteString +revcompl' = B.reverse . B.map (fromIntegral . ord . compl . chr . fromIntegral) + +compl :: Char -> Char +compl 'A' = 'T' +compl 'T' = 'A' +compl 'C' = 'G' +compl 'G' = 'C' +compl 'a' = 't' +compl 't' = 'a' +compl 'c' = 'g' +compl 'g' = 'c' +compl x = x + +defragSeq = id diff --git a/src/Config.hs b/src/Config.hs --- a/src/Config.hs +++ b/src/Config.hs @@ -13,7 +13,8 @@ import EmpFile import Version -import Bio.Sequence +import Bio.Core.Sequence +import Bio.Sequence.Fasta import Generations.GenBase import Generations.GS20 @@ -46,12 +47,12 @@ , 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)" + , flowkey = def &= explicit &= C.name "key" &= help "sequence key to start each read (TCAG)" &= typ "STRING" , 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)" + , flowlength = def &= help "number of flows for each read" &= typ "INT" + , flowcycle = def &= help "sequence nucleotides in each flow cycle (TACG)" &= typ "STRING" , inputs = def &= args &= typFile - , output = def &= help "output file" + , output = def &= help "output file" &= typFile } &= program "flowsim" &= summary ("flowsim "++version) &= details hs @@ -62,16 +63,16 @@ -- , "Trimming filters (TFILT): ...." ] -mkconf :: IO (Generation,IO [Sequence Nuc],HPLprob,FilePath) +mkconf :: IO (Generation,IO [Sequence],HPLprob,FilePath) mkconf = do cf <- cmdArgs opts - let inp = map castToNuc `fmap` if null (inputs cf) then hReadFasta stdin else concat `fmap` mapM readFasta (inputs cf) + let inp = 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) when (null $ output cf) $ error "Please specify an output file with -o" hplc <- case hplinput cf of [] -> return default_hplc - file -> (mkHPLprobs . freqtable . concatMap (hpls (f_cycle gen) . seqdata)) `fmap` readFasta file + file -> (mkHPLprobs . freqtable . concatMap (hpls (f_cycle gen) . unSD . seqdata)) `fmap` readFasta file return (gen, inp, hplc, output cf) default_hplc :: HPLprob @@ -96,8 +97,9 @@ mkgen :: Config -> IO Generation 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 + setG = setFromGen gen -- set from a 454 generation + setD = setFromStr gen -- statistical distributions + setS f s = if null s then f gen else s -- set from string value setF fs | null fs = disc_filters gen | otherwise = concatMap (maybe [] id) $ map (flip lookup filterlist) fs setT fs | null fs = trim_filters gen @@ -120,9 +122,9 @@ , models = my_model -- , disc_filters = setF df -- , trim_filters = setT tf - , f_key = setG f_key fk + , f_key = setS f_key fk , f_len = if fl > 0 then 4*((fl+3)`div` 4) else f_len gen -- is this sufficient? - , f_cycle = setG f_cycle fc + , f_cycle = setS f_cycle fc } parse_models :: String -> [Distribution] diff --git a/src/Duplicator.hs b/src/Duplicator.hs --- a/src/Duplicator.hs +++ b/src/Duplicator.hs @@ -5,16 +5,16 @@ import System.IO import System.Environment (getArgs) -import Bio.Sequence +import Bio.Sequence.Fasta import Statistics main :: IO () main = do args <- getArgs - case args of [p] -> hReadFasta stdin >>= evalRandIO . dup (read p) . map castToNuc >>= hWriteFasta stdout + case args of [p] -> hReadFasta stdin >>= evalRandIO . dup (read p) >>= 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 :: RandomGen g => Double -> [Sequence] -> Rand g [Sequence] dup p (x:xs) = do r <- getRandomR (0,1) if r < p diff --git a/src/FilterSFF.hs b/src/FilterSFF.hs new file mode 100644 --- /dev/null +++ b/src/FilterSFF.hs @@ -0,0 +1,83 @@ +{-# Language DeriveDataTypeable #-} +module FilterSFF (main) where + +-- import Debug.Trace +import System.IO + +import System.Console.CmdArgs as C +import Version + +import Generations.GenBase +import Generations.GS20 +import Generations.Titanium +import Generations.Empirical + +generations :: [(String,Generation)] +generations = [("GS20", gs20), ("Titanium",titanium),("EmpTitanium",tiEmp)] + +data Opts = Opts { generation :: String + , output, inputs :: FilePath + , count :: Bool + } + deriving (Data,Typeable,Show) + +defopt :: Opts +defopt = Opts + { generation = "Titanium" &= help "454 generation" &= typ "GEN" &= C.name "G" + , output = "filtered.sff" &= typFile + , inputs = "" &= args &= typFile + , count = False &= help "count matches instead of filtering" + } &= program "filtersff" + &= summary ("filtersff "++version) + +main :: IO () +main = do + opts <- cmdArgs defopt + SFF h rs <- readSFF (inputs opts) + case lookup (generation opts) generations of + Just g -> if count opts + then apply_count (disc_filters g) rs + else do + let f = apply_filter (disc_filters g) + n <- writeSFF' (output opts) (SFF h $ f rs) + putStrLn ("Wrote "++show n++" records to "++output opts) + Nothing -> error ("Unknown generation: '"++generation opts++"'") + +apply_filter :: [DiscardFilter] -> [ReadBlock] -> [ReadBlock] +apply_filter filters = filter (\r -> (and (apply1 filters r))) + +apply_count :: [DiscardFilter] -> [ReadBlock] -> IO () +apply_count filters rbs = go (replicate (length filters) 0) rbs + where go :: [Int] -> [ReadBlock] -> IO () + go counts (r:rs) = do + let fs = apply1 filters r + if and fs + then go counts rs + else do + let c = add counts fs + trace (show c) + go c rs + go counts [] = trace (show counts++"\n") + +{- +-- This leaks (well, retains) memory. Laziness be damned. + +apply_trace :: [DiscardFilter] -> [ReadBlock] -> [ReadBlock] +apply_trace filters = go (replicate (length filters) 0) + where go counts (r:rs) = + let fs = apply1 filters r + in if and fs then r : go counts rs + else let c = add counts fs in trace ("Filtered: "++show c) $ go counts rs + go _ [] = trace "\n" [] +-} + +trace :: String -> IO () +trace msg = hPutStr stderr ("\r"++msg) + +-- apply fs = filter (\r -> not (and $ apply1 fs r)) + +add :: [Int] -> [Bool] -> [Int] +add cs = zipWith (+) cs . map (fromEnum . not) + +apply1 :: [DiscardFilter] -> ReadBlock -> [Bool] +apply1 filters r = map ($r) filters diff --git a/src/FlowSim.hs b/src/FlowSim.hs --- a/src/FlowSim.hs +++ b/src/FlowSim.hs @@ -15,7 +15,8 @@ module Main where -import Bio.Sequence +import Bio.Core.Sequence +import Bio.Sequence.Fasta import Bio.Sequence.SFF import qualified Data.ByteString.Lazy.Char8 as B @@ -53,16 +54,15 @@ whenLoud $ inform ("Wrote "++show n++" reads to '"++o++"'.") -- testing -test_s :: Sequence Nuc -test_s = Seq (fromStr "foo") (fromStr "aacacattcgtggtnagctacggaacacattcgtggtnagctacggaacacattcgtggtnagctacggaacacattcgtggtnagctacgg") Nothing +test_s :: Sequence +test_s = Seq (SeqLabel $ B.pack "foo") (SeqData $ B.pack "aacacattcgtggtnagctacggaacacattcgtggtnagctacggaacacattcgtggtnagctacggaacacattcgtggtnagctacgg") Nothing -- | Heavy lifting. Building an SFF file from the specified information -sim454 :: RandomGen g => Generation -> HPLprob -> [Sequence Nuc] -> Rand g SFF +sim454 :: RandomGen g => Generation -> HPLprob -> [Sequence] -> Rand g SFF sim454 gen hplc ss = do let ch = makeCommonHeader gen 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 + rbs <- map tf `fmap` mapM (makeReadBlock gen hplc ch) ss return (SFF ch rbs) -- Generate a sequence of models @@ -93,15 +93,17 @@ -- generate the "absolute" flow values from the cycle and the origin sequence data 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 notElem (toUpper x) "ACGT" then B.tail s1 else s1; _ -> s1 - in c1 ++ makeFlows c s2 +makeFlows c (SeqData s) + | B.null s = [] + | otherwise = let (c1,SeqData s1) = makeCycle [] (take 4 c) (SeqData s) + s2 = case B.uncons s1 of Just (x,_) -> if notElem (toUpper x) "ACGT" then B.tail s1 else s1; _ -> s1 + in c1 ++ makeFlows c (SeqData s2) makeCycle :: [(Char,Int)] -> [Char] -> SeqData -> ([(Char,Int)],SeqData) makeCycle acc [] s = (reverse acc, s) -makeCycle acc (c:cs) s = let (this,rest) = B.span ((==toUpper c).toUpper) s - in makeCycle ((c,fromIntegral $ B.length this):acc) cs rest +makeCycle acc (c:cs) (SeqData s) = + let (this,rest) = B.span ((==toUpper c).toUpper) s + in makeCycle ((c,fromIntegral $ B.length this):acc) cs (SeqData rest) -- | Takes a sequence of flows with exact homopolymer lengths, -- permutes and calls them. Next step is prepareData @@ -148,10 +150,10 @@ -- | Generate a ReadBlock -- direction and position chosen at random, and encoded in the read name -makeReadBlock :: RandomGen g => Generation -> HPLprob -> CommonHeader -> Sequence Nuc -> Rand g ReadBlock +makeReadBlock :: RandomGen g => Generation -> HPLprob -> CommonHeader -> Sequence -> Rand g ReadBlock makeReadBlock g hplc ch sq = do let sdata = seqdata sq - rn = BS.concat $ BL.toChunks $ seqlabel sq + rn = BS.concat $ BL.toChunks $ unSL $ seqlabel sq fs = makeFlows (BS.unpack $ flow ch) sdata (pfs,cs,qs,is) <- convertCalls `fmap` permuteAndCall g hplc fs return $ verifyRB (flow_length ch) ReadBlock { @@ -167,8 +169,8 @@ -- The data block , flow_data = packFlows pfs -- :: [Flow] , flow_index = BS2.pack $ map fromIntegral is -- :: ByteString - , bases = fromStr cs -- :: SeqData - , quality = BL.pack qs -- :: QualData + , bases = SeqData $ B.pack cs + , quality = QualData $ BL.pack $ map unQual qs } -- | Consistency check on generated ReadBlocks. @@ -180,8 +182,10 @@ | 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") + | (B.length $ unSD $ bases rb) /= fromIntegral (BS.length $ flow_index rb) + = err ("bases ("++show (B.length $ unSD $ bases rb)++") and flow_index ("++show (BS.length $ flow_index rb)++") have different lengths") + | (B.length $ unQD $ quality rb) /= fromIntegral (BS.length $ flow_index rb) + = err ("quality ("++show (B.length $ unQD $ 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 --- a/src/GelFilter.hs +++ b/src/GelFilter.hs @@ -5,10 +5,11 @@ import System.Environment (getArgs) import System.IO -import Bio.Sequence +import Bio.Core.Sequence +import Bio.Sequence.Fasta main :: IO () main = do [min,max] <- map read `fmap` getArgs - hWriteFasta stdout =<< filter (\x -> seqlength x>=min && seqlength x<= max) `fmap` hReadFasta stdin + hWriteFasta stdout =<< filter (\x -> seqlength x>=(Offset min) && seqlength x<= (Offset 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 @@ -2,6 +2,7 @@ import Generations.GenBase +import Bio.Core.Sequence import Data.Array import qualified Data.ByteString.Char8 as BC import qualified Data.ByteString as B @@ -29,7 +30,7 @@ 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 = [discard_empty] + , disc_filters = [discard_key "TCAG", discard_empty] , trim_filters = [trim_sigint, trim_qual20 10] , f_key = "TCAG" @@ -40,7 +41,7 @@ qual_gs20_tab :: QualMethod -qual_gs20_tab _ _ (_,fi) = B.unpack (gs20qt!fi) +qual_gs20_tab _ _ (_,fi) = map Qual $ B.unpack (gs20qt!fi) gs20qt :: QualTable gs20qt = mkQualTable (0,3100) $ map (\(x,y) -> (x,BC.pack y)) [ diff --git a/src/Generations/GenBase.hs b/src/Generations/GenBase.hs --- a/src/Generations/GenBase.hs +++ b/src/Generations/GenBase.hs @@ -8,12 +8,17 @@ ) where import qualified Data.ByteString.Char8 as BC + +import Bio.Core.Sequence import Bio.Sequence.SFF import Bio.Sequence.SFF_filters import Statistics import HPLCount hiding (main) +-- workaround for missing Integral instance for Qual +-- round' (Qual q) = Qual (round q) + type Model = Int -> Distribution -- ^ each hpl length has its own distribution of flow values type QualMethod = HPLprob -> Model -> (Char,Flow) -> [Qual] @@ -41,6 +46,10 @@ -- ,"flow distributions:"]) -- ++ unlines [" "++show i++": "++show (model g i) | i <- [0..6]] +-- | Set absolute limits for quality values +clamp_qual :: Double -> Double -> Double +clamp_qual x = max 0 . min x + -- | Calculate exact qualities from model using Bayes' theorem -- This should emulate Marguiles et al.'s method. qual_exact_decreasing :: HPLprob -> Model -> (Char,Flow) -> [Qual] @@ -49,7 +58,7 @@ 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 = replicate (h-20) 60 ++ (map (round . max 0 . min 60 . ((-10)*) . logBase 10) $ probs) + quals = replicate (h-20) 60 ++ (map (Qual . round . clamp_qual 60 . ((-10)*) . logBase 10) $ probs) in -- trace ("#"++show (c,flow) ++"\n"++show probs ++"\n"++show quals) $ take (fromIntegral $ (fl+50) `div` 100) $ quals @@ -59,8 +68,8 @@ 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)) + prob = ft c hpl * pdf (m hpl) f / sum [pdf (m x) f * ft c x | x <- [0..20]] - 0.0001 :: Double + in replicate hpl (Qual $ round $ clamp_qual 40 $ (-10) * logBase 10 (1-prob)) -- -------------------------------------------------- -- Putting it together diff --git a/src/Generations/Titanium.hs b/src/Generations/Titanium.hs --- a/src/Generations/Titanium.hs +++ b/src/Generations/Titanium.hs @@ -10,7 +10,7 @@ , 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 = [discard_empty,discard_dots 0.05,discard_mixed,discard_length 186] -- discard_key + , disc_filters = [discard_key "TCAG", discard_empty,discard_dots 0.05,discard_mixed,discard_length 186] , trim_filters = [trim_sigint, trim_qual20 10] -- find_primer , f_key = "TCAG" , f_len = 800 -- :: Int16 diff --git a/src/HPLCount.hs b/src/HPLCount.hs --- a/src/HPLCount.hs +++ b/src/HPLCount.hs @@ -1,7 +1,8 @@ -- count homopolymer distributions in a fasta file module HPLCount where -import Bio.Sequence +import Bio.Core.Sequence +import Bio.Sequence.Fasta import qualified Data.ByteString.Lazy.Char8 as B import qualified Data.Map as M import Text.Printf @@ -20,7 +21,7 @@ main :: IO () main = do [f] <- getArgs - display "tacg" . freqtable . concatMap (hpls "tacg" . seqdata) =<< readFasta f + display "tacg" . freqtable . concatMap (hpls "tacg" . unSD . seqdata) =<< readFasta f display :: [Char] -> M.Map (Char,Int) Int -> IO () display fs m = let imax = maximum $ map snd $ M.keys m @@ -36,7 +37,7 @@ m' = M.insert i v' m in v' `seq` m' `seq` go m' is -hpls :: [Char] -> SeqData -> [(Char,Int)] +hpls :: [Char] -> B.ByteString -> [(Char,Int)] hpls xs = go (cycle xs') where xs' = map toUpper xs go (f:fs) sd | B.null sd = [] diff --git a/src/KitSim.hs b/src/KitSim.hs --- a/src/KitSim.hs +++ b/src/KitSim.hs @@ -9,21 +9,21 @@ import System.IO import qualified Data.ByteString.Lazy.Char8 as B -import Bio.Sequence +import Bio.Core.Sequence +import Bio.Sequence.Fasta main :: IO () main = do cf <- cmdArgs conf - let reader = fmap (map castToNuc) $ - case input cf of "-" -> hReadFasta stdin + let reader = 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 :: Conf -> [Sequence] -> [Sequence] apply c = map apply1 - where apply1 (Seq s d _) = Seq s (B.concat [k,d,a]) Nothing + where apply1 (Seq s (SeqData d) _) = Seq s (SeqData $ B.concat [k,d,a]) Nothing k = B.pack (key c) a = B.pack (adapter c) diff --git a/src/Mutator.hs b/src/Mutator.hs --- a/src/Mutator.hs +++ b/src/Mutator.hs @@ -7,7 +7,9 @@ import System.IO import System.Console.CmdArgs -import Bio.Sequence +import Bio.Core.Sequence +import Bio.Sequence.Fasta +import Data.ByteString.Lazy.Char8 as B import Statistics import Version @@ -36,17 +38,17 @@ x -> readFasta x outp = case output c of "-" -> hWriteFasta stdout x -> writeFasta x - inp >>= doMutate c . map castToNuc >>= outp + inp >>= doMutate c >>= outp -doMutate :: Conf -> [Sequence Nuc] -> IO [Sequence Nuc] +doMutate :: Conf -> [Sequence] -> IO [Sequence] doMutate cf = evalRandIO . mutate (subst cf) (indel cf) mutate :: RandomGen g => - Double -> Double -> [Sequence Nuc] -> Rand g [Sequence Nuc] + Double -> Double -> [Sequence] -> Rand g [Sequence] mutate sub ind = mapM mut1 where mut1 (Seq h d _) = do - d2 <- go (toStr d) - return $ Seq h (fromStr d2) Nothing + d2 <- go (B.unpack $ unSD d) + return $ Seq h (SeqData $ B.pack d2) Nothing nuc :: Int -> Char nuc x = case x of {0 -> 'a'; 1 -> 'c'; 2 -> 'g'; 3 -> 't'} go "" = return "" diff --git a/src/Statistics.hs b/src/Statistics.hs --- a/src/Statistics.hs +++ b/src/Statistics.hs @@ -8,6 +8,8 @@ ) where import Control.Monad.Random +import Control.Applicative +import Data.Char import Data.Array.Unboxed foreign import ccall "erf" erf :: Double -> Double @@ -25,11 +27,26 @@ | Empirical Double Double (UArray Int Double) -- ^ Stores a cumulative distribution, centered on first -- param (ref. 'worsen'), and second param is step size - deriving (Show,Read) + deriving (Show) -- add more: StudentsT, SkewNormal, SkewT --- dummy to make deriving Read work -instance Read (UArray i d) where +parseDist :: String -> Either String Distribution +parseDist str = case words (map toLower str) of + ("uniform":rs) -> case rs of [x,y] -> Uniform <$> num x <*> num y + _ -> Left ("Incorrect number of arguments to 'Uniform': "++show rs) + ("lognormal":rs) -> case rs of [x,y] -> LogNormal <$> num x <*> num y + _ -> Left ("Incorrect number of arguments to 'LogNormal': "++show rs) + ("normal":rs) -> case rs of [x,y] -> Normal <$> num x <*> num y + _ -> Left ("Incorrect number of arguments to 'Normal': "++show rs) + _ -> Left ("Failed to parse as distribution: '"++str++"'.") + +num :: String -> Either String Double +num str = case reads str of [(x,"")] -> Right x + _ -> Left ("Couldn't parse '"++str++"' as a number") + +instance Read Distribution where + readsPrec _ str = case parseDist str of Right d -> [(d,"")] + Left e -> error ("Couldn't understand the distribution you specified:\n"++e) -- | Build an empirical probablility distribution by mapping the given probabilities -- to uniformly spaced points starting at 'start' with 'step' points per unit. diff --git a/src/Version.hs b/src/Version.hs --- a/src/Version.hs +++ b/src/Version.hs @@ -1,4 +1,4 @@ module Version where version :: String -version = " v0.3, copyright 2010 Ketil Malde" +version = " v0.3.3, copyright 2010-2012 Ketil Malde"