packages feed

flowsim 0.3 → 0.3.3

raw patch · 16 files changed

+233/−87 lines, 16 filesdep +biocoredep +biofastadep +biosffdep −bionew-component:exe:filtersff

Dependencies added: biocore, biofasta, biosff, mtl

Dependencies removed: bio

Files

flowsim.cabal view
@@ -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 <ketil@malde.org>@@ -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
src/CloneSim.hs view
@@ -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
src/Config.hs view
@@ -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]
src/Duplicator.hs view
@@ -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
+ src/FilterSFF.hs view
@@ -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
src/FlowSim.hs view
@@ -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)
src/GelFilter.hs view
@@ -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
src/Generations/Empirical.hs view

file too large to diff

src/Generations/GS20.hs view
@@ -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)) [
src/Generations/GenBase.hs view
@@ -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
src/Generations/Titanium.hs view
@@ -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
src/HPLCount.hs view
@@ -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 = []
src/KitSim.hs view
@@ -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)
src/Mutator.hs view
@@ -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 ""
src/Statistics.hs view
@@ -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.
src/Version.hs view
@@ -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"