rbr 0.8.5 → 0.8.6
raw patch · 6 files changed
+353/−4 lines, 6 filesdep ~basedep ~bio
Dependency ranges changed: base, bio
Files
- rbr.cabal +5/−4
- src/Lib/FreqTable.lhs +94/−0
- src/Lib/Options.lhs +76/−0
- src/Lib/Stats.hs +136/−0
- src/Lib/TextFormat.lhs +24/−0
- src/Lib/Util.hs +18/−0
rbr.cabal view
@@ -1,5 +1,5 @@ Name: rbr-Version: 0.8.5+Version: 0.8.6 License: GPL License-File: LICENSE Stability: Beta@@ -13,11 +13,10 @@ . The Darcs repository is at <http://malde.org/~ketil/biohaskell/rbr>. Homepage: http://malde.org/~ketil/- Cabal-Version: >= 1.2 Tested-With: GHC==6.10.4 -- add fps for ghc-6.4.2-Build-Depends: base>3 && <4.2, bio >= 0.4, QuickCheck, binary, bytestring, containers+Build-Depends: base>3 && <5, bio >= 0.4, QuickCheck, binary, bytestring, containers Build-Type: Simple data-files: idea.txt, README, TODO@@ -25,14 +24,16 @@ Executable rbr Executable: rbr Hs-Source-Dirs: src+ Other-Modules: Lib.Options, Lib.Stats, Lib.Util, Lib.FreqTable Main-Is: RBR.lhs- build-depends: base >3&&<4.2, bytestring, containers, bio >= 0.4+ build-depends: base >3 && <5, bytestring, containers, bio >= 0.3.1 Extensions: CPP Ghc-Options: -Wall -funbox-strict-fields Executable mct Executable: mct Hs-Source-Dirs: src+ Other-Modules: Lib.TextFormat Main-Is: MCT.lhs build-depends: bytestring Extensions: BangPatterns, FlexibleContexts
+ src/Lib/FreqTable.lhs view
@@ -0,0 +1,94 @@+\section{FreqTable}++Building a Frequency Table for keys.++Constructs the table of all possible keys, so be careful.++\begin{code}++module Lib.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 a] -> 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 a] -> 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 a] -> 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 a] -> 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}
+ src/Lib/Options.lhs view
@@ -0,0 +1,76 @@+\subsection{Option handling}++You know the drill...++\begin{code}++module Lib.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.6"+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 ['M'] ["mask-sample"] (ReqArg (\arg opt -> return opt { masksample = Just arg }) "file") "Sample to use for building masking index"+ ,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 = arg : lib opt}) "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 :: [FilePath]+ , skip :: Int+ , masksample :: Maybe FilePath+ }++defaultopts :: Opts+defaultopts = Opts 16 Nothing True False 0 0 0 Lower False False False [] 0 Nothing++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}
+ src/Lib/Stats.hs view
@@ -0,0 +1,136 @@+{-# LANGUAGE CPP, Rank2Types #-}+module Lib.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 '*')++-}
+ src/Lib/TextFormat.lhs view
@@ -0,0 +1,24 @@+\begin{code}++-- todo: support uneven column lengths++module Lib.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}
+ src/Lib/Util.hs view
@@ -0,0 +1,18 @@+module Lib.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