SFF_filters.hs

-- | This implements a number of filters used in the Titanium pipeline
module Bio.Sequence.SFF_filters where

import Bio.Sequence.SFF (ReadBlock(..), ReadHeader(..), flowToBasePos, flowgram)
import qualified Data.ByteString.Lazy as B
import qualified Data.ByteString.Lazy.Char8 as BL
import Data.List (tails)

-- Ti uses a set of filters, described in the (something) manual.
-- (GS Run Processor Application, section 3.2.2++)

-- ** Discarding filters **

-- | DiscardFilters determine whether a read is to be retained or discarded
type DiscardFilter = ReadBlock -> Bool -- True to retain, False to discard

filter_dots, filter_mixed, filter_key, filter_empty :: DiscardFilter

filter_empty rb = num_bases (read_header rb) >= 5

filter_key rb = ("TCAG"==) $ take 4 $ BL.unpack $ bases rb

-- | 3.2.2.1.2 The "dots" filter discards sequences where the last positive flow is 
--   before flow 84, and flows with >5% dots (i.e. three successive noise values) 
--   before the last postitive flow. (Interpreted as 5% of called sequence length is Ns?)
filter_dots rb = let dots = BL.length $ BL.filter (=='N') $ bases rb
                 in fromIntegral dots / fromIntegral (BL.length $ bases rb) < (0.05 :: Double)

-- | 3.2.2.1.3 The "mixed" filter discards sequences with more than 70% positive flows.  
--   Also, discard with <30% noise, >20% middle (0.45..0.75) or <30% positive.
filter_mixed rb = let fs = dropWhile (<50) . reverse . flowgram $ rb
                      fl = dlength fs
                  in and [ (dlength (filter (>50) fs) / fl) > 0.7
                         , (dlength (filter (<45) fs) / fl) > 0.3 
                         , (dlength (filter (>75) fs) / fl) > 0.3
                         , (dlength (filter (\f -> f<=75 && f>=45) fs) / fl) < 0.2
                         ]
-- | Discard a read if the number of untrimmed flows is less than n (n=186 for Titanium)
filter_length :: Int -> DiscardFilter
filter_length n rb = length (flowgram rb) >= n

-- ** Trimming filters **

-- | TrimFilters modify the read, typically trimming it for quality
type TrimFilter = ReadBlock -> ReadBlock
filter_sigint, filter_qual20 :: TrimFilter

-- | 3.2.2.1.4 Signal intensity trim - trim back until <3% borderline flows (0.5..0.7).
--   Then trim borderline values or dots from the end (use a window).
filter_sigint rb = clipFlows rb (sigint rb)

sigint :: ReadBlock -> Int
sigint rb = let bs = scanl (\(n,m,_) f -> if f >= 50 && f <= 70 then (n+1,m+1,f) else (n,m+1,f)) (0,0,0) $ flowgram rb 
            in length $ reverse 
                   $ dropWhile (\(_,_,f) -> f<=70)
                   $ dropWhile (\(n,m,_)->n `div` (m*1000) > (30::Int)) $ reverse bs


-- 3.2.2.1.5 Primer filter and 3.2.2.1.6 Trimback valley filter are ignored, 
-- since we don't know how to detect the primer, and don't understand the description of tbv.

-- | 3.2.2.1.7 Quality score trimming trims using a 10-base window until a Q20 average is found. 
filter_qual20 rs = clipSeq rs $ qual20 rs

qual20 :: ReadBlock -> Int
qual20 rs = (fromIntegral $ num_bases $ read_header rs) 
            - (length . takeWhile (<20) . map (avg . take 10) . tails . reverse . B.unpack $ quality rs)

-- ** Utility functions **

-- | List length as a double (eliminates many instances of fromIntegral)
dlength :: [a] -> Double
dlength = fromIntegral . length

-- | Calculate average of a list
avg :: Integral a => [a] -> Double
avg xs = fromIntegral (sum xs) / dlength xs

-- | Translate a number of flows to position in sequence, and update clipping data accordingly
clipFlows :: ReadBlock -> Int -> ReadBlock
clipFlows rb n = clipSeq rb (flowToBasePos rb n)

-- | Update clip_qual_right if more severe than previous value
clipSeq :: ReadBlock -> Int -> ReadBlock
clipSeq rb n' = let n = fromIntegral n' 
                    rh = read_header rb 
                in if clip_qual_right rh <= n then rb else rb { read_header = rh {clip_qual_right = n }}