module Generations.GenBase ( Model, QualMethod , Generation(..), makeCommonHeader , qual_exact_fixed, qual_exact_decreasing , module Bio.Sequence.SFF , module Bio.Sequence.SFF_filters , module Statistics , module HPLCount ) where import qualified Data.ByteString.Char8 as BC import Bio.Sequence.SFF import Bio.Sequence.SFF_filters import Debug.Trace (trace) import Statistics import HPLCount hiding (main) type Model = Int -> Distribution -- ^ each hpl length has its own distribution of flow values type QualMethod = HPLprob -> Model -> (Char,Flow) -> [Qual] data Generation = Gen { name :: String , qcall :: QualMethod -- ^ use the current model and hpl lenght to generate qual data , models :: [Model] -- ^ a progression of initial models for the sequence of flows , degrade :: Distribution -- the amount to increase a stdev multiplier each position , disc_filters :: [DiscardFilter] , trim_filters :: [TrimFilter] -- stuff to build the commonheader from , f_key :: String , f_len :: Int , f_cycle :: String -- readblock :: [ReadBlock] ? } instance Show Generation where show g = "Config:\n" ++ unlines (map (" "++) ["key: "++f_key g ,"flowlength: "++show (f_len g), "flow cycle: "++show (f_cycle g) ]) -- ,"flow distributions:"]) -- ++ unlines [" "++show i++": "++show (model g i) | i <- [0..6]] -- | 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]] probs = drop 1 $ scanl (+) 0 $ map prob [0..20] quals = 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 -- | 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 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)) -- -------------------------------------------------- -- Putting it together makeCommonHeader :: Generation -> CommonHeader makeCommonHeader g = CommonHeader { index_offset = 0 , index_length = 0 , key_length = fromIntegral $ length $ f_key g , flowgram_fmt = 1 , key = BC.pack $ f_key g , num_reads = 0 -- fromIntegral $ n_reads g , flow_length = fromIntegral $ f_len g , flow = BC.pack $ concat $ replicate (f_len g `div` (length $ f_cycle g)) $ f_cycle g }