GenBase.hs

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.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]

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] ?
    , adapter :: String
    }

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]]

-- | 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]
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 = 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
     
-- | 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,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]] - 0.0001 :: Double
  in replicate hpl (Qual $ round $ clamp_qual 40 $ (-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
    }