rbr 0.8.3 → 0.8.5
raw patch · 12 files changed
+505/−868 lines, 12 filesdep ~basedep ~bionew-uploader
Dependency ranges changed: base, bio
Files
- RBR/FreqTable.lhs +0/−94
- RBR/MCT.lhs +0/−105
- RBR/Options.lhs +0/−74
- RBR/RBR.lhs +0/−385
- RBR/Stats.hs +0/−136
- RBR/TextFormat.lhs +0/−24
- RBR/Util.hs +0/−18
- RBR/hooks.c +0/−15
- README +2/−2
- rbr.cabal +10/−15
- src/MCT.lhs +105/−0
- src/RBR.lhs +388/−0
− RBR/FreqTable.lhs
@@ -1,94 +0,0 @@-\section{FreqTable}--Building a Frequency Table for keys.--Constructs the table of all possible keys, so be careful.--\begin{code}--module RBR.FreqTable (FreqTable(..)- ,freqtable_int,freqtable_integer -- constructors- ,sparsetable_int, sparsetable_integer- ) where--import Bio.Sequence--import Data.IntMap as IM-import Data.Map as M-import Data.List as L--data FreqTable a = FT { count :: a -> Int- , counts :: [(a,Int)]- }--\end{code}--\subsection{Constructing the frequency table}--Frequency table, counting occurences.--\begin{code}--freqtable_int :: HashF Int -> [Sequence] -> FreqTable Int-freqtable_int kf ss = FT (\k -> IM.findWithDefault 0 k idx) (IM.toList idx)- where- idx = foldl' myupdate IM.empty $ concatMap (L.map fst . hashes kf . seqdata) ss- myupdate m k = let x = 1 + IM.findWithDefault 0 k m- in x `seq` IM.insert k x m--freqtable_integer :: HashF Integer -> [Sequence] -> FreqTable Integer-freqtable_integer kf ss = FT (\k -> M.findWithDefault 0 k idx) (M.toList idx)- where- idx = foldl' myupdate M.empty $ concatMap (L.map fst . hashes kf . seqdata) ss- myupdate m k = let x = 1 + M.findWithDefault 0 k m- in x `seq` M.insert k x m----- alternative: use--- accumArray (+) 0 (minBound, maxBound) . (map (\x->(x,1)))--\end{code}--\subsection{Sparse maps}--For parameter k, increment all existing keys, insert new keys only-when a gap larger than k would otherwise result. Keys are added-weight according to distance to next inserted key, so that the sum-weight of a sequence in the map is independent of k.--For this to be useful, masking must be performed against a sliding-average. (How will this work against libraries?)--Invariant: sparsetable_int 1 == freqtable_int--\begin{code}--sparsetable_int :: Int -> HashF Int -> [Sequence] -> FreqTable Int-sparsetable_int step kf ss = FT (\k -> find k idx) (IM.toList idx)- where- find = IM.findWithDefault 0- idx = foldl' ins1 IM.empty ss- ins1 ix s = let ks = L.map fst $ hashes kf $ seqdata s- in go ix (step-1) ks- -- insert keys that already exist, or when the skip-counter is zero- go ix' n (k:ks) | n==0 || find k ix' /= 0- || L.null ks = let x = find k ix'+step-n- in x `seq` go (IM.insert k x ix') (step-1) ks- | True = go ix' (n-1) ks- go ix' _n [] = ix'--sparsetable_integer :: Int -> HashF Integer -> [Sequence] -> FreqTable Integer-sparsetable_integer step kf ss = FT (\k -> find k idx) (M.toList idx)- where- find = M.findWithDefault 0- idx = foldl' ins1 M.empty ss- ins1 ix s = let ks = L.map fst $ hashes kf $ seqdata s- in go ix (step-1) ks- -- insert keys that already exist, or when the skip-counter is zero- go ix' n (k:ks) | n==0 || find k ix' /= 0- || L.null ks = let x = find k ix'+step-n- in x `seq` go (M.insert k x ix') (step-1) ks- | True = go ix' (n-1) ks- go ix' _n [] = ix'--\end{code}
− RBR/MCT.lhs
@@ -1,105 +0,0 @@-\section{MaskCount}--Examine each nucleotide from two data sets, and count the Ns (or Xs?).-Separate into A AB B 0. Also (optionally?) count masked region sizes-to generate gnuplot'able histograms.--TODO:- - parameter to specify lower case/N/X as masked- - justify columns in output- - output percentages--\begin{code}-{-# LANGUAGE BangPatterns, FlexibleContexts #-}-module Main where--import Bio.Sequence--import RBR.TextFormat--import Control.Monad (when)-import Data.Char (isLower)-import Data.List (foldl',intersperse,partition)-import System.Environment (getArgs)-import System.Exit (exitWith,ExitCode(..))-import Text.Printf (printf, PrintfType)--main :: IO ()-main = do- (opts,[a,b]) <- parseArgs- as <- readFasta a- bs <- readFasta b- let res = zipComp as bs- when (elem 'l' opts)- (do- putStrLn ("# fraction of masked positions: "++a++" both "++b++"none")- putStrLn $ unlines $ map printres res)- when (elem 's' opts)- (do- putStrLn "# Summary:"- putStrLn $ printSum a b $ foldl' add (CS 0 0 0 0) (map snd res))--parseArgs :: IO (String, [String])-parseArgs = do- as <- getArgs- let (opts,files) = partition ((=='-').head) as- when (length files /= 2) usage- let short = concat $ map tail opts- return (short,files)--usage :: IO a-usage = do putStrLn "mc: compare two FASTA files for masked nucleotides"- putStrLn "usage: mc [-s|-l] file_1 file_2"- exitWith $ ExitFailure 1--data CompStruct = CS !Int !Int !Int !Int--add :: CompStruct -> CompStruct -> CompStruct-add (CS a ab b o) (CS a' ab' b' o') = CS (a+a') (ab+ab') (b+b') (o+o')---- construct the percentage similar nucleotides-printres :: (String, CompStruct) -> String-printres (s,CS a ab b o) = let tot = a+ab+b+o- in concat $ intersperse " "- $ s : map (percent tot) [a,ab,b,o]--printSum :: String -> String -> CompStruct -> String-printSum a b (CS ca cab cb co) = unlines (["Summary for "++a++" vs. "++b++":\n"]- ++ columns " " [ map (" "++) ["unmasked in both:","masked in both:"- ,"masked in "++a++":","masked in "++b++":"]- , map show [co, cab, ca, cb]- , map (paren . percent tot) [co,cab,ca,cb ]])- ++ "Pearson's phi: "++printf "%.3f" (phi :: Double)- where tot = co+cab+ca+cb- phi = let [a,b,ab,o] = map ((/(fromIntegral tot)) . fromIntegral) [ca,cb,cab,co]- in (ab*o-a*b)/sqrt((ab+a)*(ab+b)*(a+o)*(b+o))--paren :: String -> String-paren s = "("++s++"%)"--percent :: (PrintfType (Double -> t), Integral a, Integral a1) => a1 -> a -> t-percent tot x = printf "%.2f" (100*fromIntegral x/fromIntegral tot :: Double)--zipComp :: [Sequence] -> [Sequence] -> [(String,CompStruct)]-zipComp as bs = map (z1 (CS 0 0 0 0)) $ zipSeqs as bs--z1 :: CompStruct -> (Sequence, Sequence) -> (String, CompStruct)-z1 cs (a,b) = (toStr $ seqlabel a, foldl' collect cs $ zip (get a) (get b))- where get s = map isLower $ map (s!) [0..seqlength s-1]- -- s/isLower/isN/ for N detect.--zipSeqs :: [Sequence] -> [Sequence] -> [(Sequence,Sequence)]-zipSeqs [] [] = []-zipSeqs (a:as) (b:bs) =- if seqlabel a == seqlabel b then (a,b):zipSeqs as bs- else error ("sequence sets differ (offending sequences: '"- ++ toStr (seqlabel a) ++ "' and '" ++ toStr (seqlabel b) ++"'.")-zipSeqs _ _ = error "extraneous sequences"--collect :: CompStruct -> (Bool, Bool) -> CompStruct-collect (CS a ab b o) (True,True) = (CS a (ab+1) b o)-collect (CS a ab b o) (True,False) = (CS (a+1) ab b o)-collect (CS a ab b o) (False,True) = (CS a ab (b+1) o)-collect (CS a ab b o) (False,False) = (CS a ab b (o+1))--\end{code}
− RBR/Options.lhs
@@ -1,74 +0,0 @@-\subsection{Option handling}--You know the drill...--\begin{code}--module RBR.Options (version,usage,getOptions,Opts(..),parseargs,MaskWith(..)) where--import System.Environment (getArgs)-import System.Console.GetOpt-import System.Exit--version, usagemsg :: String-version = "0.8.4"-usagemsg = "Usage: rbr [-k <k>] [options] FILE.."--getOptions :: IO ([Opts -> IO Opts], [String], [String])-getOptions = getArgs >>= (return . getOpt Permute options)--usage :: [String] -> String-usage errs = usageInfo (concat errs ++ usagemsg) options--data MaskWith = Lower | Ns | Xs--options :: [OptDescr (Opts -> IO Opts)]-options = [Option ['k'] ["word-size"] (ReqArg (\arg opt -> return opt { kval = read arg }) "Int") "Word size"- ,Option [] ["sparse-keys"] (ReqArg (\arg opt -> return opt { skip = read arg }) "Int") "Skip length"--- ,Option [] ["shape"] (ReqArg ...) "Key shape"- ,Option ['N'] ["no-stats"] (NoArg (\opt -> return opt { stats = False })) "Omit statistical masking"- ,Option ['g'] ["gapclosing"] (ReqArg (\arg opt -> return opt { gap = read arg }) "Int") "Close gaps"- ,Option ['s'] ["stringency"] (ReqArg (\arg opt -> return opt { string = read arg }) "Float") "Mask stringecy"- ,Option ['t'] ["threshold"] (ReqArg (\arg opt -> return opt { thresh = read arg }) "Float") "Mask threshold"- ,Option [] ["preserve-lower-case"] (NoArg (\opt -> return opt { preserve_lower = True })) "Preserve lower case in input sequences"-- ,Option ['D'] ["distribution"] (NoArg (\opt -> return opt { distr = True })) "Distribution output"- ,Option ['l'] ["library"] (ReqArg (\arg opt -> return opt {lib = Just arg}) "file") "Mask words from library"- ,Option ['L'] ["lower-case"] (NoArg (\opt -> return opt { lower = Lower })) "Mask using lower case letters"- ,Option ['n'] ["n-mask"] (NoArg (\opt -> return opt { lower = Ns })) "Mask using 'n's"- ,Option ['x'] ["x-mask"] (NoArg (\opt -> return opt { lower = Xs })) "Mask using 'X's"- ,Option [] ["server"] (NoArg (\opt -> return opt { server = True })) "Run in server mode"- ,Option ['o'] ["output-file"] (ReqArg (\arg opt -> return opt { ofile = Just arg }) "file") "Output file"- ,Option ['v'] ["verbose"] (NoArg (\opt -> return opt { verb = True })) "Display progress"-- ,Option ['h'] ["help"] (NoArg (\_ -> do {putStrLn $ usage []; exitWith ExitSuccess })) "Display help"- ,Option ['V'] ["version"] (NoArg (\_ -> do {putStrLn $ "rbr version "++version; exitWith ExitSuccess })) "Display version"- ]--data Opts = Opts { kval :: Int, ofile :: Maybe FilePath- , stats, preserve_lower :: Bool- , thresh, string :: Double- , gap :: Int- , lower :: MaskWith- , distr, verb, server :: Bool- , lib :: Maybe FilePath- , skip :: Int- }--defaultopts :: Opts-defaultopts = Opts 16 Nothing True False 0 0 0 Lower False False False Nothing 0--parseargs :: [Opts -> IO Opts] -> IO Opts-parseargs args = foldl (>>=) (return defaultopts) args >>= fixS >>= fixT- where fixS opts = return $ case (string opts,lower opts) of -- sensible(?) defaults- (0,Ns) -> opts { string = 2.0 }- (0,Xs) -> opts { string = 2.0 }- (0,Lower) -> opts { string = 1.1 }- _ -> opts- fixT opts = return $ case (thresh opts,lower opts) of -- sensible(?) defaults- (0,Ns) -> opts { thresh = 8.0 }- (0,Xs) -> opts { thresh = 8.0 }- (0,Lower) -> opts { thresh = 5.0 }- _ -> opts--\end{code}
− RBR/RBR.lhs
@@ -1,385 +0,0 @@-\section{RBR - RepeatBeater revisited}--RBR works by examining word frequencies in EST data, and-based on the assumption that expression of positions within each gene-is relatively uniform, identifies extremely common words as potential-repeats.--A predecessor to RBR was implemented as some scripts parsing-\texttt{xsact} output, this is a revised, standalone implmentation.--\begin{code}-{-# LANGUAGE CPP #-}--module Main (main) where--import Bio.Sequence--import RBR.FreqTable-import RBR.Stats hiding (median,uv_add)-import RBR.Util (interactl)-import RBR.Options---- import Text.Printf-import qualified Data.ByteString.Lazy.Char8 as B (map)--import System.IO-import System.IO.Unsafe-import Control.Monad-import Data.List as List-import Data.Char as Char--import Data.Map hiding (null,filter,map,elems,keys)--#define perror (\s->error ("Program error in function \""++s++"\" at "++__FILE__++":"++show (__LINE__::Int)++\-"\nPlease report this incident to <ketil@ii.uib.no>."))--main :: IO ()-main = do- (opt1,non,err) <- getOptions- -- print (opt1,non,err)- when ((null non && null opt1) || not (null err)) (error $ usage err)- opts <- parseargs opt1- case non of [] -> do ss <- hReadFasta stdin- main_real opts ("<stdin>", return ss)- [f1] -> main_real opts (f1, readFasta f1)- fs -> error ("Only supply one input file (you supplied: "++show fs++")")--header :: Opts -> [FilePath] -> String-header opts _inputs = "# rbr version "++version++", mask="++(case (lower opts) of Lower -> "L"; Ns -> "n"; Xs -> "X")++- " k="++show (kval opts)++" stringency="++show (string opts)++- " thresh="++show (thresh opts)++" g="++show (gap opts) ++ " skip="++ show (skip opts)--\end{code}--\subsection{The real main}--\begin{code}--debug :: Bool -> String -> IO ()-debug p m = when p $ hPutStr stderr m--countIO :: String -> String -> Int -> [a] -> IO [a]-countIO msg post step xs = sequence $ map unsafeInterleaveIO ((blank >> outmsg (0::Int) >> c):cs)- where (c:cs) = ct 0 xs- output = hPutStr stderr- blank = output ('\r':take 70 (repeat ' '))- outmsg x = output ('\r':msg++show x) >> hFlush stderr- ct s ys = let (a,b) = splitAt (step-1) ys- next = s+step- in case b of [b1] -> map return a ++ [outmsg (s+step) >> hPutStr stderr post >> return b1]- [] -> map return (init a) ++ [outmsg (s+length a) >> hPutStr stderr post >> return (last a)]- _ -> map return a ++ [outmsg s >> return (head b)] ++ ct next (tail b)--main_real :: Opts -> (String, IO [Sequence]) -> IO ()-main_real opts (ifile,ssio) = do- when (verb opts || distr opts) $ putStrLn $ header opts [ifile]- when (verb opts) $ do- debug (verb opts) ("Reading file: "++ifile++"\n")- ss <- ssio- debug (verb opts) ("Number of sequences: "++show (length ss)++"\n")- if kval opts <= 16 then main_cont make_ft_int opts ssio- else main_cont make_ft_integer opts ssio--main_cont :: Integral a => MkFT a -> Opts -> IO [Sequence] -> IO ()-main_cont make_ft opts ssio = do- mask1 <- mask_func make_ft opts ssio- ss'' <- if (verb opts) then countIO " ...masking sequences: " ", done.\n" 100 =<< ssio- else ssio- if (server opts) then do- ft <- make_ft "Building index array" (verb opts) (skip opts) (kval opts) ssio- interactl "READY: " (concatMap (main_server opts ft ss'') . words)- else do if (distr opts) then main_coverage make_ft opts ssio- else main_masked opts mask1 ss''--type MkFT a = String -> Bool -> Int -> Int -> IO [Sequence] -> IO (FreqTable a)--make_ft_int :: MkFT Int-make_ft_int msg vb skp k ssio = do- ss' <- if vb then countIO (" "++msg++": ") ", done.\n" 100 =<< ssio- else ssio- if skp > 0 then return (sparsetable_int skp (rcontig k) ss')- else return (freqtable_int (rcontig k) ss')--make_ft_integer :: MkFT Integer-make_ft_integer msg vb skp k ssio = do- ss' <- if vb then countIO (" "++msg++": ") ", done.\n" 100 =<< ssio- else ssio- if skp>0 then return (sparsetable_integer skp (rcontig k) ss')- else return (freqtable_integer (rcontig k) ss')--main_coverage :: Integral a => MkFT a -> Opts -> IO [Sequence] -> IO ()-main_coverage make_ft opts ssio = do- ft <- make_ft "Building index array" (verb opts) (skip opts) (kval opts) ssio- let k = kval opts- sk = skip opts- (stdevs,b2,strg) = (thresh opts, 5, string opts)- seq_cov s = let cv = coverage k sk (rcontig k) ft s- ls = sortv $ blunt b2- $ map (\l->(fromIntegral $ head l,length l)) $ group $ sort cv- (mu,stdv) = distrib strg ls- t = mu + max b2 (stdevs*stdv)- in (seqheader s, if null cv || null ls then (0,0,0) else (mu,stdv,t),cv) -- todo: round!- write = case ofile opts of Nothing -> putStrLn; Just f -> writeFile f- debug (verb opts) "Calculating coverage\n"- ss <- ssio- write $ unlines $ map show $ map seq_cov ss--main_masked :: Opts -> (Sequence->Sequence) -> [Sequence] -> IO ()-main_masked opts mf ss = do- let ms = if kval opts == 0 then ss else map mf ss- write = case ofile opts of Nothing -> hWriteFasta stdout- Just f -> writeFasta f- write ms--mask_func :: Integral a => MkFT a -> Opts -> IO [Sequence] -> IO (Sequence -> Sequence)-mask_func make_ft opts ssio =- do let k = kval opts- g = gap opts- sk = skip opts- gc = if g > 0 then closegaps (g+k-1) else id- sparm = (thresh opts, 5, string opts)- upcase = if preserve_lower opts then id else upcaseSequence- ms <- do case stats opts of- False -> return []- True -> do ft <- make_ft "Indexing input" (verb opts) sk k ssio- return [mask_stat k sk sparm ft]- mt <- case lib opts of- Nothing -> return []- Just f -> do lft <- make_ft "Indexing repeat library" (verb opts) sk k $ readFasta f- return [mask_table k sk lft]- return ((case ms++mt of [] -> id- xs -> mask_generic (lower opts) (k+if sk==0 then 0 else sk-1)- (gc . min_mask ((sk+2) `div` 2) . combine xs)) . upcase)--upcaseSequence :: Sequence -> Sequence-upcaseSequence (Seq l d q) = (Seq l (B.map toUpper d) q)---- sparse keys can give aberrant counts, so require more than one consequtive--- todo: really really convert to RLL-encoded whitelist-min_mask :: Int -> [Bool] -> [Bool]-min_mask _ [] = []-min_mask sk xs@(True:_) = let (t,f) = span id xs in t ++ min_mask sk f-min_mask sk xs@(False:_) = let (f,t) = span not xs in (if length f >= sk then f- else take (length f) (repeat True)) ++ min_mask sk t--combine :: [a -> [Bool]] -> (a -> [Bool])-combine [] = perror "combine"-combine [f] = f-combine (f1:f2:fs) = let f = \s -> zipWith (&&) (f1 s) (f2 s) in combine (f:fs)---main_server :: Integral a => Opts -> FreqTable a -> [Sequence] -> (String -> String)-main_server opts ft ss =- let sidx = fromList $ map (\s->(seqlabel s,s)) ss- (k,sk,gc,ps) = (kval opts, skip opts,if gap opts > 0 then closegaps (k+gap opts-1) else id, (thresh opts,5,string opts))- in (\s -> case Data.Map.lookup (fromStr s) sidx of- Nothing -> ("input error: "++s++" not found\n")- Just sq -> let cover = show $ coverage k sk (rcontig k) ft sq- masked_ = show $ mask_generic (lower opts) (if sk==0 then k else k+sk-1) (gc . mask_stat k (skip opts) ps ft) sq- -- mask_stat (lower opts) k g ps ft sq- unmasked = show $ seqdata sq- in unlines [masked_, cover,unmasked])--\end{code}--Generic masking function. The list of Bool is a whitelist, so that-True corresponds to retained nucleotides, while False indicates-nucleotides to be masked.--(TODO: Run-length encode the boolean list?)--\begin{code}--mask_generic :: MaskWith -> Int -> (Sequence -> [Bool]) -> Sequence -> Sequence-mask_generic b sz fn s@(Seq l sd mq) = Seq l (fromStr . masked b sz (fn s) . toStr $ sd) mq--\end{code}--Mask all words from a dictionary (i.e. a library of known repeats)--\begin{code}--mask_table :: Integral a => Int -> Int -> FreqTable a -> Sequence -> [Bool]-mask_table k skp rlib s = let- cv = coverage k skp (rcontig k) rlib s -- slight abuse of 'coverage'- in map (==0) cv--\end{code}--Masking sequences from frequency counts.--\begin{code}--type MPars = (Double,Double,Double) -- stdevs thres, 1-elim, stringency---- calculate the 'whitelist'-mask_stat :: Integral a => Int -> Int -> MPars -> FreqTable a -> Sequence -> [Bool]-mask_stat k sk (stdevs,b2,strg) ft s =- let cv = coverage k sk (rcontig k) ft s- ls = sortv $ blunt b2- $ map (\l->(fromIntegral $ head l,length l)) $ group $ sort cv- (mu,stdv) = distrib strg ls- t = mu + max b2 (stdevs*stdv)- in if null cv || null ls then take (1-k+fromIntegral (seqlength s)) $ repeat True- else map ((<=t) . fromIntegral) cv---- distrib calculates the base distribution-distrib :: Double -> [(Double,Int)] -> (Double,Double)-distrib strg ls = distr' ls (uv_mk [])- where- distr' [] _- = perror "distrib"- distr' ((mag,cnt):ls') ustd =- let ustd' = uv_add mag cnt ustd- u = uniVar ustd'- score = variance u / stdev u- in if score >= strg || null ls'- then (mean u,stdev u)- else distr' ls' ustd'---- s {seqannot = seqannot s--- ++[UKV "Threshold" $ printf "%.2f" t]--- ,seqdata = listArray (bounds $ seqdata s) ms}---- compensate for 1% error rate by reducing number of ones by k%--- this is inaccurate, but possibly a workable estimate nonetheless-blunt :: Double -> [(Double,Int)] -> [(Double,Int)]-blunt k = reduce1s . dropWhile ((==0).fst)- where reduce1s all_@((a,x):ls) = if a==1- then (1,max 0 (x-floor k*(sum $ map snd all_) `div` 100)):ls- else all_- reduce1s [] = [] -- sequence contains no words---- sort ls around the centre or the modal interval-sortv :: [(Double,Int)] -> [(Double,Int)]-sortv ls = let m = scaled_modal ls- in sortBy (\(a,_) (c,_) -> compare (abs (a-m)) (abs (c-m))) ls---- NB: 1 = error - chose mode ignoring 1s? (see C.e. seq. 24)---- for each window of occurrence frequencies (n..sqrt n), find the one with most positions-scaled_modal :: [(Double,Int)] -> Double-scaled_modal = modal . map window . tails'--window :: [(Double,Int)] -> (Double,Int)-window [] = perror "window"-window ((mag,cnt):ls) =- let top = mag + sqrt mag- xs = (mag,cnt) : takeWhile ((<top).fst) ls- in ((mag+(fst . last) xs)/2, sum (map snd xs))---- tails without the final [[]]-tails' :: [a] -> [[a]]-tails' [] = []-tails' x@(_:xs) = x : tails' xs--{--windowize ((mag,count):ls) =- let (this,rest) = break ((> mag+sqrt mag) . fst) ls- in (mag+sqrt mag,count+sum (map snd this)) : windowize rest-windowize [] = []--}--modal :: Ord b => [(a,b)] -> a-modal = fst . maximumBy (\x y -> compare (snd x) (snd y))--median :: [(Double,Int)] -> Double-median ls = let n = (`div` 2) $ sum $ map snd ls- in nquant n ls--nquant :: Int -> [(Double,Int)] -> Double-nquant n ((v,c):vs) | n < c = v- | otherwise = nquant (n-c) vs-nquant _ [] = perror "median of an empty list is impossible!"----- faster(?) replacements for Stats.uv and uv_del-type UV = (Int,Double,Double,Double,Double)--uv_mk :: [(Double,Int)] -> UV-uv_mk = foldr uv1 (0,0,0,0,0)- where uv1 (m,cnt) (n,x,x2,x3,x4) = let c = fromIntegral cnt in- (n+cnt,x+m*c,x2+m*m*c,x3+m*m*m*c,x4+m*m*m*m*c)--uv_sub, uv_add :: Double -> Int -> UV -> UV-uv_sub m cnt (n,x,x2,x3,x4) = let c = fromIntegral cnt in- (n-cnt,x-m*c,x2-m*m*c,x3-m*m*m*c,x4-m*m*m*m*c)--uv_add m cnt (n,x,x2,x3,x4) = let c = fromIntegral cnt in- (n+cnt,x+m*c,x2+m*m*c,x3+m*m*m*c,x4+m*m*m*m*c)---- masked takes mask type, kval, threshold, a whitelist ('False' means mask)--- and the (nucleotides in the) sequence-masked :: MaskWith -> Int -> [Bool] -> String -> String-masked Lower = masked' toLower 0-masked Ns = masked' (const 'n') 0-masked Xs = masked' (const 'X') 0---- mask using an arbitrary masking function--- the 'ns' parameter keeps track of any remaining masked word, so that--- we finish masking out the whole over-represented word-masked' :: (Char->Char) -> Int -> Int -> [Bool] -> String -> String-masked' _ _ _ (_:_) [] =- perror "masked' (1)"-masked' mf ns _ [] es = map mf (take ns es) ++ drop ns es-masked' mf ns k (b:bs) (e:es)- | ns == 0 && b = e : masked' mf 0 k bs es- | ns > 0 && b = mf e : masked' mf (ns-1) k bs es- | not b = mf e : masked' mf k k bs es- | otherwise = perror "masked' (2)"---- close gaps of less than 'g' unmasked words between groups of masked characters--- 'False' corresponds to words to mask (whitelist)-closegaps :: Int -> [Bool] -> [Bool]-closegaps g = close g True . partitions---- always starts with a (possibly empty) unmasked group-partitions :: [Bool] -> [Int]-partitions es = map List.length $ if head es == False then []:partitions' es else partitions' es- where partitions' = List.groupBy (\x y -> x == False && y == False || x /= False && y /= False)---- convert runs of <=g Trues to Falses-close :: Int -> Bool -> [Int] -> [Bool]-close _ _ [] = []-close g True (c:cs) = take c (repeat (if c>g then True else False)) ++ close g False cs-close g False (c:cs) = (take c $ repeat False) ++ close g True cs--\end{code}--\subsection{Coverage}--Calculating the coverage over a sequence--For sparse indices, this should be a sliding average.--\begin{code}--coverage :: Integral a => Int -> Int -> HashF a -> FreqTable a -> Sequence -> [Int]-coverage k skp kf tb s = if skp > 0 then sliding_avg skp raw_counts else raw_counts- where- raw_counts = map mcount $ padkeys (seqlength s-fromIntegral k+1) 0 $ hashes kf $ seqdata s- mcount (Just k1) = count tb k1- mcount Nothing = 0--padkeys :: Offset -> Offset -> [(k,Offset)] -> [Maybe k]-padkeys top cur [] = take (fromIntegral (top - cur)) $ repeat Nothing-padkeys top cur ((k,p):ks) = if cur == p then Just k : padkeys top (cur+1) ks- else Nothing : padkeys top (cur+1) ((k,p):ks)----- or just:--- sliding_avg k = map ((`div` k) . sum . take k) . tails--sliding_avg :: Int -> [Int] -> [Int]-sliding_avg s xs = let x1 = take s xs- in if length x1 == s then slav s (sum x1) xs (drop s xs)- else [sum x1 `div` length x1]---- todo: use floats?-slav :: Int -> Int -> [Int] -> [Int] -> [Int]-slav s tmp hd tl@(_:_) = tmp `div` s : slav s (tmp-head hd+head tl) (tail hd) (tail tl)-slav s tmp _ [] = [tmp `div` s]--\end{code}-
− RBR/Stats.hs
@@ -1,136 +0,0 @@-{-# LANGUAGE CPP, Rank2Types #-}-module RBR.Stats (UniVar(..),uniVar,uv,uv_add,uv_del- ,Quantiles(..),quantiles) where -- ,histogram,display) where--import Data.List (sort, sortBy, group)--- import FiniteMap--class Statistic s where- samples :: s -> Int--data UniVar = UV { uv_samples :: Int- , mean :: Double- , stdev :: Double- , variance :: Double- , skewness :: Double- , kurtosis :: Double- , sumSq :: Double- , coeffVar :: Double- , stdErrMn :: Double- }--instance Statistic UniVar where samples = uv_samples--instance Show UniVar where- show u = adjust- [["Samples", (show $ samples u)]- , ["Mean", (show $ mean u)]- , ["Standard dev", (show $ stdev u)]- , ["Variance", (show $ variance u)]- , ["Skewness", (show $ skewness u)]- , ["Kurtosis", (show $ kurtosis u)]- , ["Sum of squares", (show $ sumSq u)]- , ["Coeff. of var", (show $ coeffVar u)]- , ["Std err mean", (show $ stdErrMn u)]- ]--adjust :: [[String]] -> String-adjust [] = []-adjust ([a,b]:xs) = (a++":"++take (15-length a) (repeat ' ')++b) ++ "\n" ++ adjust xs-adjust _ = error "Pattern matching error in 'adjust'"--type UVTMP = (Int,Double,Double,Double,Double)---- | more or less the univariate function from SAS--- calculate by tracking n, sum of x, of x², x^3, x^4-uniVar :: (Int, Double, Double, Double, Double) -> UniVar-uniVar (n',x,x2,x3,x4) =- let n = fromIntegral n'- m = x/n- m2 = m*m- m3 = m*m2- var = (x2-m2*n)/(n-1)- s = sqrt(var)- skew = (x3 - 3*m*x2 + 2*m3*n)/(s*s*s*n)- kurt = (x4 - 4*m*x3 + 6*m2*x2 - 4*m3*x + n*m*m3)/(s*s*s*s*n) - 3- cv = s/m- in UV n' m s var skew kurt x2 cv (s/sqrt n)--uv_add, uv_del :: Real a => a -> UVTMP -> UVTMP-uv_add = uv_upd (+)-uv_del = uv_upd (-)---- requires -fglasgow-exts-uv_upd :: (Real a) => (forall b . Real b => b -> b -> b)- -> a -> UVTMP -> UVTMP-uv_upd f d' (n,x,x2,x3,x4) = let d = toFloat d' in n `seq` x `seq` x2- `seq` x3 `seq` x4 `seq` (f n 1,f x d,f x2 (d*d),f x3 (d*d*d),f x4 (d*d*d*d))- where toFloat = fromRational . toRational---uv :: Real a => [a] -> UVTMP-uv ds = foldr uv_add (0,0,0,0,0) ds--data Quantiles = Qs { wsamples :: Int- , smallest :: Double- , quartile1 :: Double- , median :: Double- , mode :: [Double]- , quartile3 :: Double- , greatest :: Double- }--instance Statistic Quantiles where samples = wsamples--instance Show Quantiles where- show w = adjust- [["Samples", (show $ samples w)]- , ["Smallest", (show $ smallest w)]- , ["Q1", (show $ quartile1 w)]- , ["Median", (show $ median w)]- , ["Modes", (show $ mode w)]- , ["Q3", (show $ quartile3 w)]- , ["Greatest", (show $ greatest w)]]--quantiles :: [Double] -> Quantiles-quantiles ds = let- n = length ds- sorted = sort ds- q1 = case n `quotRem` 4 of (q,0) -> ((sorted!!(q-1))+(sorted!!q))/2.0- (q,_) -> sorted!!q- q2 = case n `quotRem` 2 of (q,0) -> ((sorted!!(q-1))+(sorted!!q))/2.0- (q,_) -> sorted!!q- q3 = case (3*n) `quotRem` 4 of (q,0) -> ((sorted!!(q-1))+(sorted!!q))/2.0- (q,_) -> sorted!!q- modes = let- ms = sortOn (negate.fst) $ map (\x->(length x,head x))$ group sorted- in (snd $ head ms) : map snd- (takeWhile (\x->fst x==fst (head ms)) (tail ms))- in- Qs n (head sorted) q1 q2 modes q3 (last sorted)--sortOn :: Ord b => (a->b) -> [a] -> [a]-sortOn f = sortBy (\x y -> compare (f x) (f y))--{--type Histogram = FiniteMap Double Int---- | histogram builds a histogram given the list of midpoints-histogram :: [Double] -> [Double] -> Histogram-histogram ms xs = foldl (insert ms) emptyFM' xs- where- emptyFM' = foldl (\s v -> addToFM s v 0) emptyFM ms- insert :: [Double] -> Histogram -> Double -> Histogram- insert (m1:m2:ms) s x = if abs (m1-x) <= abs (m2-x)- then addToFM_C (+) s m1 1 else insert (m2:ms) s x- insert [m1] s x = addToFM_C (+) s m1 1- insert [] _ _ = error "Must provide at least one midpoint"---- todo: speed up with strict foldl'--display :: Histogram -> String-display h = unlines $ map disp1 $ fmToList h- where disp1 (v,n) = show v ++ (take (7-(length $ show v)) (repeat ' ')) ++- ": " ++ (take n $ repeat '*')---}
− RBR/TextFormat.lhs
@@ -1,24 +0,0 @@-\begin{code}---- todo: support uneven column lengths--module RBR.TextFormat where--import Data.List--columns :: String -> [[String]] -> [String]-columns delim = map (concat . intersperse delim) . transpose . extend--extend :: [[String]] -> [[String]]-extend cs = zipWith pad (sizes cs) cs---- mkrow rsizes = concat . intersperse "^" . zipWith pad rsizes--pad :: Int -> [String] -> [String]-pad n = map (pad1 n)- where pad1 o text = text ++ take (o - length text) (repeat ' ')--sizes :: [[String]] -> [Int]-sizes = map (maximum . map length)--\end{code}
− RBR/Util.hs
@@ -1,18 +0,0 @@-module RBR.Util where--import System.IO--interactl :: String -> (String -> String) -> IO ()-interactl p f = do- putStr p- hFlush stdout- ls <- getContents- i' p f (lines ls)--i' :: String -> (t -> String) -> [t] -> IO ()-i' _ _ [] = return ()-i' p f (l:ls) = do- putStr (f l)- putStr p- hFlush stdout- i' p f ls
− RBR/hooks.c
@@ -1,15 +0,0 @@-/* after a tip from David Roundy */-#include <Rts.h>-#include <RtsFlags.h>-#include <unistd.h>--void defaultsHook (void) {- RtsFlags.GcFlags.maxStkSize = 8*10000002 / sizeof(W_); /* 80M */--#ifdef _SC_PHYS_PAGES- unsigned long long pagesize = sysconf(_SC_PAGESIZE);- unsigned long long numpages = sysconf(_SC_PHYS_PAGES);- unsigned long long mhs = numpages*pagesize*8/10;- RtsFlags.GcFlags.maxHeapSize = 1ULL+mhs/BLOCK_SIZE_W;-#endif-}
README view
@@ -30,7 +30,7 @@ currently works. You need GHC (http://haskell.org/ghc). Everything is tested against-version 6.6, but version 6.4.x might also work. Expect to do some+version 6.8, but older versions might also work. Expect to do some modifications in the code for earlier versions than that. You first need to get my 'bio' library, it is available from the same@@ -40,7 +40,7 @@ chmod +x Setup.hs ./Setup.hs configure (add --prefix=$HOME if you don't have root access) ./Setup.hs build- ./Setup.hs install+ {sudo} ./Setup.hs install If you want to go the more manual route, cd to the src subdirectory, and
rbr.cabal view
@@ -1,5 +1,5 @@ Name: rbr-Version: 0.8.3+Version: 0.8.5 License: GPL License-File: LICENSE Stability: Beta@@ -11,34 +11,29 @@ statistical model to identify suspicious sequence parts, and masks them. The README has more details. .- The Darcs repository is at <http://malde.org/~ketil/rbr>.+ The Darcs repository is at <http://malde.org/~ketil/biohaskell/rbr>. Homepage: http://malde.org/~ketil/ Cabal-Version: >= 1.2-Tested-With: GHC==6.8.2+Tested-With: GHC==6.10.4 -- add fps for ghc-6.4.2-Build-Depends: base, bio >= 0.3.1, QuickCheck, binary, bytestring, containers+Build-Depends: base>3 && <4.2, bio >= 0.4, QuickCheck, binary, bytestring, containers Build-Type: Simple data-files: idea.txt, README, TODO -Library- Exposed-Modules: RBR.Options, RBR.Util, RBR.Stats, RBR.FreqTable, RBR.TextFormat- C-sources: RBR/hooks.c- Build-Depends: base>3, containers, bio- Extensions:- Ghc-Options: -Wall -funbox-strict-fields- Executable rbr Executable: rbr- Main-Is: RBR/RBR.lhs- build-depends: bytestring+ Hs-Source-Dirs: src+ Main-Is: RBR.lhs+ build-depends: base >3&&<4.2, bytestring, containers, bio >= 0.4 Extensions: CPP Ghc-Options: -Wall -funbox-strict-fields Executable mct Executable: mct- Main-Is: RBR/MCT.lhs+ Hs-Source-Dirs: src+ Main-Is: MCT.lhs build-depends: bytestring- Extensions: BangPatterns, FlexibleContexts+ Extensions: BangPatterns, FlexibleContexts Ghc-Options: -Wall -funbox-strict-fields
+ src/MCT.lhs view
@@ -0,0 +1,105 @@+\section{MaskCount}++Examine each nucleotide from two data sets, and count the Ns (or Xs?).+Separate into A AB B 0. Also (optionally?) count masked region sizes+to generate gnuplot'able histograms.++TODO:+ - parameter to specify lower case/N/X as masked+ - justify columns in output+ - output percentages++\begin{code}+{-# LANGUAGE BangPatterns, FlexibleContexts #-}+module Main where++import Bio.Sequence++import Lib.TextFormat++import Control.Monad (when)+import Data.Char (isLower)+import Data.List (foldl',intersperse,partition)+import System.Environment (getArgs)+import System.Exit (exitWith,ExitCode(..))+import Text.Printf (printf, PrintfType)++main :: IO ()+main = do+ (opts,[a,b]) <- parseArgs+ as <- readFasta a+ bs <- readFasta b+ let res = zipComp as bs+ when (elem 'l' opts)+ (do+ putStrLn ("# fraction of masked positions: "++a++" both "++b++"none")+ putStrLn $ unlines $ map printres res)+ when (elem 's' opts)+ (do+ putStrLn "# Summary:"+ putStrLn $ printSum a b $ foldl' add (CS 0 0 0 0) (map snd res))++parseArgs :: IO (String, [String])+parseArgs = do+ as <- getArgs+ let (opts,files) = partition ((=='-').head) as+ when (length files /= 2) usage+ let short = concat $ map tail opts+ return (short,files)++usage :: IO a+usage = do putStrLn "mc: compare two FASTA files for masked nucleotides"+ putStrLn "usage: mc [-s|-l] file_1 file_2"+ exitWith $ ExitFailure 1++data CompStruct = CS !Int !Int !Int !Int++add :: CompStruct -> CompStruct -> CompStruct+add (CS a ab b o) (CS a' ab' b' o') = CS (a+a') (ab+ab') (b+b') (o+o')++-- construct the percentage similar nucleotides+printres :: (String, CompStruct) -> String+printres (s,CS a ab b o) = let tot = a+ab+b+o+ in concat $ intersperse " "+ $ s : map (percent tot) [a,ab,b,o]++printSum :: String -> String -> CompStruct -> String+printSum a b (CS ca cab cb co) = unlines (["Summary for "++a++" vs. "++b++":\n"]+ ++ columns " " [ map (" "++) ["unmasked in both:","masked in both:"+ ,"masked in "++a++":","masked in "++b++":"]+ , map show [co, cab, ca, cb]+ , map (paren . percent tot) [co,cab,ca,cb ]])+ ++ "Pearson's phi: "++printf "%.3f" (phi :: Double)+ where tot = co+cab+ca+cb+ phi = let [a,b,ab,o] = map ((/(fromIntegral tot)) . fromIntegral) [ca,cb,cab,co]+ in (ab*o-a*b)/sqrt((ab+a)*(ab+b)*(a+o)*(b+o))++paren :: String -> String+paren s = "("++s++"%)"++percent :: (PrintfType (Double -> t), Integral a, Integral a1) => a1 -> a -> t+percent tot x = printf "%.2f" (100*fromIntegral x/fromIntegral tot :: Double)++zipComp :: [Sequence a] -> [Sequence a] -> [(String,CompStruct)]+zipComp as bs = map (z1 (CS 0 0 0 0)) $ zipSeqs as bs++z1 :: CompStruct -> (Sequence a, Sequence a) -> (String, CompStruct)+z1 cs (a,b) = (toStr $ seqlabel a, foldl' collect cs $ zip (get a) (get b))+ where get s = map isLower $ map (s!) [0..seqlength s-1]+ -- s/isLower/isN/ for N detect.++zipSeqs :: [Sequence a] -> [Sequence a] -> [(Sequence a,Sequence a)]+zipSeqs [] [] = []+zipSeqs (a:as) (b:bs) =+ if seqlabel a == seqlabel b then (a,b):zipSeqs as bs+ else error ("sequence sets differ (offending sequences: '"+ ++ toStr (seqlabel a) ++ "' and '" ++ toStr (seqlabel b) ++"'.")+zipSeqs _ _ = error "extraneous sequences"++collect :: CompStruct -> (Bool, Bool) -> CompStruct+collect (CS a ab b o) (True,True) = (CS a (ab+1) b o)+collect (CS a ab b o) (True,False) = (CS (a+1) ab b o)+collect (CS a ab b o) (False,True) = (CS a ab (b+1) o)+collect (CS a ab b o) (False,False) = (CS a ab b (o+1))++\end{code}
+ src/RBR.lhs view
@@ -0,0 +1,388 @@+\section{RBR - RepeatBeater revisited}++RBR works by examining word frequencies in EST data, and+based on the assumption that expression of positions within each gene+is relatively uniform, identifies extremely common words as potential+repeats.++A predecessor to RBR was implemented as some scripts parsing+\texttt{xsact} output, this is a revised, standalone implmentation.++\begin{code}+{-# LANGUAGE CPP #-}++module Main (main) where++import Bio.Sequence++import Lib.FreqTable+import Lib.Stats hiding (median,uv_add)+import Lib.Util (interactl)+import Lib.Options++-- import Text.Printf+import qualified Data.ByteString.Lazy.Char8 as B (map)++import System.IO+import System.IO.Unsafe+import Control.Monad+import Data.List as List+import Data.Char as Char++import Data.Map hiding (null,filter,map,elems,keys)++#define perror (\s->error ("Program error in function \""++s++"\" at "++__FILE__++":"++show (__LINE__::Int)++\+"\nPlease report this incident to <ketil@ii.uib.no>."))++main :: IO ()+main = do+ (opt1,non,err) <- getOptions+ -- print (opt1,non,err)+ when ((null non && null opt1) || not (null err)) (error $ usage err)+ opts <- parseargs opt1+ case non of [] -> do ss <- return . map castToNuc =<< hReadFasta stdin+ main_real opts ("<stdin>", return ss)+ [f1] -> main_real opts (f1, return . map castToNuc =<< readFasta f1)+ fs -> error ("Only supply one input file (you supplied: "++show fs++")")++header :: Opts -> [FilePath] -> String+header opts _inputs = "# rbr version "++version++", mask="++(case (lower opts) of Lower -> "L"; Ns -> "n"; Xs -> "X")+++ " k="++show (kval opts)++" stringency="++show (string opts)+++ " thresh="++show (thresh opts)++" g="++show (gap opts) ++ " skip="++ show (skip opts)++\end{code}++\subsection{The real main}++\begin{code}++debug :: Bool -> String -> IO ()+debug p m = when p $ hPutStr stderr m++countIO :: String -> String -> Int -> [a] -> IO [a]+countIO msg post step xs = sequence $ map unsafeInterleaveIO ((blank >> outmsg (0::Int) >> c):cs)+ where (c:cs) = ct 0 xs+ output = hPutStr stderr+ blank = output ('\r':take 70 (repeat ' '))+ outmsg x = output ('\r':msg++show x) >> hFlush stderr+ ct s ys = let (a,b) = splitAt (step-1) ys+ next = s+step+ in case b of [b1] -> map return a ++ [outmsg (s+step) >> hPutStr stderr post >> return b1]+ [] -> map return (init a) ++ [outmsg (s+length a) >> hPutStr stderr post >> return (last a)]+ _ -> map return a ++ [outmsg s >> return (head b)] ++ ct next (tail b)++main_real :: Opts -> (String, IO [Sequence Nuc]) -> IO ()+main_real opts (ifile,ssio) = do+ when (verb opts || distr opts) $ putStrLn $ header opts [ifile]+ when (verb opts) $ do+ debug (verb opts) ("Reading file: "++ifile++"\n")+ ss <- ssio+ debug (verb opts) ("Number of sequences: "++show (length ss)++"\n")+ if kval opts <= 16 then main_cont make_ft_int opts ssio+ else main_cont make_ft_integer opts ssio++main_cont :: Integral a => MkFT a -> Opts -> IO [Sequence Nuc] -> IO ()+main_cont make_ft opts ssio = do+ mask1 <- mask_func make_ft opts ssio+ ss'' <- if (verb opts) then countIO " ...masking sequences: " ", done.\n" 100 =<< ssio+ else ssio+ if (server opts) then do+ ft <- make_ft "Building index array" (verb opts) (skip opts) (kval opts) ssio+ interactl "READY: " (concatMap (main_server opts ft ss'') . words)+ else do if (distr opts) then main_coverage make_ft opts ssio+ else main_masked opts mask1 ss''++type MkFT a = String -> Bool -> Int -> Int -> IO [Sequence Nuc] -> IO (FreqTable a)++make_ft_int :: MkFT Int+make_ft_int msg vb skp k ssio = do+ ss' <- if vb then countIO (" "++msg++": ") ", done.\n" 100 =<< ssio+ else ssio+ if skp > 0 then return (sparsetable_int skp (rcontig k) ss')+ else return (freqtable_int (rcontig k) ss')++make_ft_integer :: MkFT Integer+make_ft_integer msg vb skp k ssio = do+ ss' <- if vb then countIO (" "++msg++": ") ", done.\n" 100 =<< ssio+ else ssio+ if skp>0 then return (sparsetable_integer skp (rcontig k) ss')+ else return (freqtable_integer (rcontig k) ss')++main_coverage :: Integral a => MkFT a -> Opts -> IO [Sequence Nuc] -> IO ()+main_coverage make_ft opts ssio = do+ ft <- make_ft "Building index array" (verb opts) (skip opts) (kval opts) ssio+ let k = kval opts+ sk = skip opts+ (stdevs,b2,strg) = (thresh opts, 5, string opts)+ seq_cov s = let cv = coverage k sk (rcontig k) ft s+ ls = sortv $ blunt b2+ $ map (\l->(fromIntegral $ head l,length l)) $ group $ sort cv+ (mu,stdv) = distrib strg ls+ t = mu + max b2 (stdevs*stdv)+ in (seqheader s, if null cv || null ls then (0,0,0) else (mu,stdv,t),cv) -- todo: round!+ write = case ofile opts of Nothing -> putStrLn; Just f -> writeFile f+ debug (verb opts) "Calculating coverage\n"+ ss <- ssio+ write $ unlines $ map show $ map seq_cov ss++main_masked :: Opts -> (Sequence a->Sequence a) -> [Sequence a] -> IO ()+main_masked opts mf ss = do+ let ms = if kval opts == 0 then ss else map mf ss+ write = case ofile opts of Nothing -> hWriteFasta stdout+ Just f -> writeFasta f+ write ms++mask_func :: Integral a => MkFT a -> Opts -> IO [Sequence Nuc] -> IO (Sequence Nuc -> Sequence Nuc)+mask_func make_ft opts ssio =+ do let k = kval opts+ g = gap opts+ sk = skip opts+ gc = if g > 0 then closegaps (g+k-1) else id+ sparm = (thresh opts, 5, string opts)+ upcase = if preserve_lower opts then id else upcaseSequence+ ms <- do case masksample opts of + Nothing -> case stats opts of+ False -> return []+ True -> do ft <- make_ft "Indexing input" (verb opts) sk k ssio+ return [mask_stat k sk sparm ft]+ Just f -> do ft <- make_ft "Indexing sample" (verb opts) sk k (return . map castToNuc =<< readFasta f)+ return [mask_stat k sk sparm ft]+ mt <- mapM (\libf -> do+ lft <- make_ft "Indexing repeat library" (verb opts) sk k (return . map castToNuc =<< readFasta libf)+ return $ mask_table k sk lft) (lib opts) + return ((case ms++mt of [] -> id+ xs -> mask_generic (lower opts) (k+if sk==0 then 0 else sk-1)+ (gc . min_mask ((sk+2) `div` 2) . combine xs)) . upcase)++upcaseSequence :: Sequence a -> Sequence a+upcaseSequence (Seq l d q) = (Seq l (B.map toUpper d) q)++-- sparse keys can give aberrant counts, so require more than one consequtive+-- todo: really really convert to RLL-encoded whitelist+min_mask :: Int -> [Bool] -> [Bool]+min_mask _ [] = []+min_mask sk xs@(True:_) = let (t,f) = span id xs in t ++ min_mask sk f+min_mask sk xs@(False:_) = let (f,t) = span not xs in (if length f >= sk then f+ else take (length f) (repeat True)) ++ min_mask sk t++combine :: [a -> [Bool]] -> (a -> [Bool])+combine [] = perror "combine"+combine [f] = f+combine (f1:f2:fs) = let f = \s -> zipWith (&&) (f1 s) (f2 s) in combine (f:fs)+++main_server :: Integral a => Opts -> FreqTable a -> [Sequence Nuc] -> (String -> String)+main_server opts ft ss =+ let sidx = fromList $ map (\s->(seqlabel s,s)) ss+ (k,sk,gc,ps) = (kval opts, skip opts,if gap opts > 0 then closegaps (k+gap opts-1) else id, (thresh opts,5,string opts))+ in (\s -> case Data.Map.lookup (fromStr s) sidx of+ Nothing -> ("input error: "++s++" not found\n")+ Just sq -> let cover = show $ coverage k sk (rcontig k) ft sq+ masked_ = show $ mask_generic (lower opts) (if sk==0 then k else k+sk-1) (gc . mask_stat k (skip opts) ps ft) sq+ -- mask_stat (lower opts) k g ps ft sq+ unmasked = show $ seqdata sq+ in unlines [masked_, cover,unmasked])++\end{code}++Generic masking function. The list of Bool is a whitelist, so that+True corresponds to retained nucleotides, while False indicates+nucleotides to be masked.++(TODO: Run-length encode the boolean list?)++\begin{code}++mask_generic :: MaskWith -> Int -> (Sequence Nuc -> [Bool]) -> Sequence Nuc -> Sequence Nuc+mask_generic b sz fn s@(Seq l sd mq) = Seq l (fromStr . masked b sz (fn s) . toStr $ sd) mq++\end{code}++Mask all words from a dictionary (i.e. a library of known repeats)++\begin{code}++mask_table :: Integral a => Int -> Int -> FreqTable a -> Sequence Nuc -> [Bool]+mask_table k skp rlib s = let+ cv = coverage k skp (rcontig k) rlib s -- slight abuse of 'coverage'+ in map (==0) cv++\end{code}++Masking sequences from frequency counts.++\begin{code}++type MPars = (Double,Double,Double) -- stdevs thres, 1-elim, stringency++-- calculate the 'whitelist'+mask_stat :: Integral a => Int -> Int -> MPars -> FreqTable a -> Sequence Nuc -> [Bool]+mask_stat k sk (stdevs,b2,strg) ft s =+ let cv = coverage k sk (rcontig k) ft s+ ls = sortv $ blunt b2+ $ map (\l->(fromIntegral $ head l,length l)) $ group $ sort cv+ (mu,stdv) = distrib strg ls+ t = mu + max b2 (stdevs*stdv)+ in if null cv || null ls then take (1-k+fromIntegral (seqlength s)) $ repeat True+ else map ((<=t) . fromIntegral) cv++-- distrib calculates the base distribution+distrib :: Double -> [(Double,Int)] -> (Double,Double)+distrib strg ls = distr' ls (uv_mk [])+ where+ distr' [] _+ = perror "distrib"+ distr' ((mag,cnt):ls') ustd =+ let ustd' = uv_add mag cnt ustd+ u = uniVar ustd'+ score = variance u / stdev u+ in if score >= strg || null ls'+ then (mean u,stdev u)+ else distr' ls' ustd'++-- s {seqannot = seqannot s+-- ++[UKV "Threshold" $ printf "%.2f" t]+-- ,seqdata = listArray (bounds $ seqdata s) ms}++-- compensate for 1% error rate by reducing number of ones by k%+-- this is inaccurate, but possibly a workable estimate nonetheless+blunt :: Double -> [(Double,Int)] -> [(Double,Int)]+blunt k = reduce1s . dropWhile ((==0).fst)+ where reduce1s all_@((a,x):ls) = if a==1+ then (1,max 0 (x-floor k*(sum $ map snd all_) `div` 100)):ls+ else all_+ reduce1s [] = [] -- sequence contains no words++-- sort ls around the centre or the modal interval+sortv :: [(Double,Int)] -> [(Double,Int)]+sortv ls = let m = scaled_modal ls+ in sortBy (\(a,_) (c,_) -> compare (abs (a-m)) (abs (c-m))) ls++-- NB: 1 = error - chose mode ignoring 1s? (see C.e. seq. 24)++-- for each window of occurrence frequencies (n..sqrt n), find the one with most positions+scaled_modal :: [(Double,Int)] -> Double+scaled_modal = modal . map window . tails'++window :: [(Double,Int)] -> (Double,Int)+window [] = perror "window"+window ((mag,cnt):ls) =+ let top = mag + sqrt mag+ xs = (mag,cnt) : takeWhile ((<top).fst) ls+ in ((mag+(fst . last) xs)/2, sum (map snd xs))++-- tails without the final [[]]+tails' :: [a] -> [[a]]+tails' [] = []+tails' x@(_:xs) = x : tails' xs++{-+windowize ((mag,count):ls) =+ let (this,rest) = break ((> mag+sqrt mag) . fst) ls+ in (mag+sqrt mag,count+sum (map snd this)) : windowize rest+windowize [] = []+-}++modal :: Ord b => [(a,b)] -> a+modal = fst . maximumBy (\x y -> compare (snd x) (snd y))++median :: [(Double,Int)] -> Double+median ls = let n = (`div` 2) $ sum $ map snd ls+ in nquant n ls++nquant :: Int -> [(Double,Int)] -> Double+nquant n ((v,c):vs) | n < c = v+ | otherwise = nquant (n-c) vs+nquant _ [] = perror "median of an empty list is impossible!"+++-- faster(?) replacements for Stats.uv and uv_del+type UV = (Int,Double,Double,Double,Double)++uv_mk :: [(Double,Int)] -> UV+uv_mk = foldr uv1 (0,0,0,0,0)+ where uv1 (m,cnt) (n,x,x2,x3,x4) = let c = fromIntegral cnt in+ (n+cnt,x+m*c,x2+m*m*c,x3+m*m*m*c,x4+m*m*m*m*c)++uv_sub, uv_add :: Double -> Int -> UV -> UV+uv_sub m cnt (n,x,x2,x3,x4) = let c = fromIntegral cnt in+ (n-cnt,x-m*c,x2-m*m*c,x3-m*m*m*c,x4-m*m*m*m*c)++uv_add m cnt (n,x,x2,x3,x4) = let c = fromIntegral cnt in+ (n+cnt,x+m*c,x2+m*m*c,x3+m*m*m*c,x4+m*m*m*m*c)++-- masked takes mask type, kval, threshold, a whitelist ('False' means mask)+-- and the (nucleotides in the) sequence+masked :: MaskWith -> Int -> [Bool] -> String -> String+masked Lower = masked' toLower 0+masked Ns = masked' (const 'n') 0+masked Xs = masked' (const 'X') 0++-- mask using an arbitrary masking function+-- the 'ns' parameter keeps track of any remaining masked word, so that+-- we finish masking out the whole over-represented word+masked' :: (Char->Char) -> Int -> Int -> [Bool] -> String -> String+masked' _ _ _ (_:_) [] =+ perror "masked' (1)"+masked' mf ns _ [] es = map mf (take ns es) ++ drop ns es+masked' mf ns k (b:bs) (e:es)+ | ns == 0 && b = e : masked' mf 0 k bs es+ | ns > 0 && b = mf e : masked' mf (ns-1) k bs es+ | not b = mf e : masked' mf k k bs es+ | otherwise = perror "masked' (2)"++-- close gaps of less than 'g' unmasked words between groups of masked characters+-- 'False' corresponds to words to mask (whitelist)+closegaps :: Int -> [Bool] -> [Bool]+closegaps g = close g True . partitions++-- always starts with a (possibly empty) unmasked group+partitions :: [Bool] -> [Int]+partitions [] = [] -- at least mask_table may return an empty list here+partitions es = map List.length $ if head es == False then []:partitions' es else partitions' es+ where partitions' = List.groupBy (\x y -> x == False && y == False || x /= False && y /= False)++-- convert runs of <=g Trues to Falses+close :: Int -> Bool -> [Int] -> [Bool]+close _ _ [] = []+close g True (c:cs) = take c (repeat (if c>g then True else False)) ++ close g False cs+close g False (c:cs) = (take c $ repeat False) ++ close g True cs++\end{code}++\subsection{Coverage}++Calculating the coverage over a sequence++For sparse indices, this should be a sliding average.++\begin{code}++coverage :: Integral a => Int -> Int -> HashF a -> FreqTable a -> Sequence Nuc -> [Int]+coverage k skp kf tb s = if skp > 0 then sliding_avg skp raw_counts else raw_counts+ where+ raw_counts = map mcount $ padkeys (seqlength s-fromIntegral k+1) 0 $ hashes kf $ seqdata s+ mcount (Just k1) = count tb k1+ mcount Nothing = 0++padkeys :: Offset -> Offset -> [(k,Offset)] -> [Maybe k]+padkeys top cur [] = take (fromIntegral (top - cur)) $ repeat Nothing+padkeys top cur ((k,p):ks) = if cur == p then Just k : padkeys top (cur+1) ks+ else Nothing : padkeys top (cur+1) ((k,p):ks)+++-- or just:+-- sliding_avg k = map ((`div` k) . sum . take k) . tails++sliding_avg :: Int -> [Int] -> [Int]+sliding_avg s xs = let x1 = take s xs+ in if length x1 == s then slav s (sum x1) xs (drop s xs)+ else [sum x1 `div` length x1]++-- todo: use floats?+slav :: Int -> Int -> [Int] -> [Int] -> [Int]+slav s tmp hd tl@(_:_) = tmp `div` s : slav s (tmp-head hd+head tl) (tail hd) (tail tl)+slav s tmp _ [] = [tmp `div` s]++\end{code}+