diff --git a/rbr.cabal b/rbr.cabal --- a/rbr.cabal +++ b/rbr.cabal @@ -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 . 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 diff --git a/src/Lib/FreqTable.lhs b/src/Lib/FreqTable.lhs new file mode 100644 --- /dev/null +++ b/src/Lib/FreqTable.lhs @@ -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} diff --git a/src/Lib/Options.lhs b/src/Lib/Options.lhs new file mode 100644 --- /dev/null +++ b/src/Lib/Options.lhs @@ -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 ] [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} diff --git a/src/Lib/Stats.hs b/src/Lib/Stats.hs new file mode 100644 --- /dev/null +++ b/src/Lib/Stats.hs @@ -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 '*') + +-} diff --git a/src/Lib/TextFormat.lhs b/src/Lib/TextFormat.lhs new file mode 100644 --- /dev/null +++ b/src/Lib/TextFormat.lhs @@ -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} diff --git a/src/Lib/Util.hs b/src/Lib/Util.hs new file mode 100644 --- /dev/null +++ b/src/Lib/Util.hs @@ -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