clustertools 0.1.1 → 0.1.2
raw patch · 2 files changed
+72/−40 lines, 2 filesdep +QuickCheck
Dependencies added: QuickCheck
Files
- clustertools.cabal +2/−2
- src/ClusterLibs.hs +70/−38
clustertools.cabal view
@@ -1,5 +1,5 @@ Name: clustertools-Version: 0.1.1+Version: 0.1.2 License: GPL License-file: LICENSE Author: Ketil Malde@@ -72,7 +72,7 @@ Executable clusterlibs Hs-source-dirs: src Main-Is: ClusterLibs.hs- Build-depends: regex-compat+ Build-depends: regex-compat, QuickCheck Executable xcerpt Hs-source-dirs: src
src/ClusterLibs.hs view
@@ -2,34 +2,51 @@ -- calculate clusters by library, using the lib table -- read patterns from a table, first line is header -{--import qualified Data.ByteString.Char8.Lazy as B-import Data.Map (Map)-import qualified Data.Map as M---type Str = B.ByteString-type Clusters = Map Str [Str]--insert1 m k v = M.insertWith (flip (++)) k [v] m--}- import System.Environment (getArgs) import Text.Regex import Data.List (elemIndex,sortBy,intersperse)+import Data.Map hiding (map)+import Numeric (showFFloat) +import Statistics++type LibTable = [(Regex,String)] type Cluster = (String,[String]) main :: IO () main = do [ps,cs] <- do [xs,ys] <- getArgs return [xs,ys] - `catch` error "Usage: clusterlibs ps cs"+ `catch` error "Usage: clusterlibs <libtable> <clusters>" pat <- readPatternTable ps+ counts <- totalsByLib pat cs+ writeTable [("","sum":map snd pat++["significance"])]+ writeTable . map (decorate counts) . countClusters pat =<< readClusters cs++-- counts is +decorate :: (String,[Int]) -> (String,[Int]) -> (String,[String])+decorate counts (x,ys) = + let fractions = [ fromIntegral x / fromIntegral (sum $ snd counts) | x <- snd counts]+ clustersize = sum ys+ significance = unwords $ zipWith isSignificant ys fractions+ isSignificant obs frac = showFFloat (Just 3) (pvalue (1-frac) clustersize (clustersize-obs)) ""+ in (x,map show (clustersize:ys) ++ [significance])++{- clus <- classClusters pat cs- writeClusters clus+ writeTable clus+-} -readPatternTable :: FilePath -> IO [(Regex,String)]+-- --------------------------------------------------+-- | Calculate the p value of observing k sequences+-- from a library, given an a priori background distribution. +pvalue :: Double -> Int -> Int -> Double+pvalue fraction tot obs = cumbin fraction tot obs++-- --------------------------------------------------+-- | Read a LibTable from a file of wsp-separated columns, first line is+-- a header, and only columns "Pattern" and "Name" are used.+readPatternTable :: FilePath -> IO LibTable readPatternTable f = do (h:ls) <- return . map words . lines =<< readFile f let (p,n) = case (elemIndex "Pattern" h, elemIndex "Name" h) of@@ -39,34 +56,49 @@ | otherwise = (mkRegex (l!!p),l!!n) return $ map z ls --- will need to match against all, to check for multiple matches--- tag names with library-classClusters :: [(Regex,String)] -> FilePath -> IO [Cluster]-classClusters ps f = return . getClusters (classify ps) . lines =<< readFile f+-- --------------------------------------------------+-- | Acquire total counts by library from a clustering+totalsByLib :: LibTable -> FilePath -> IO (String,[Int])+totalsByLib lt f = do + cs <- readClusters f+ let unify cs = [("total", concatMap snd cs)]+ case countClusters lt (unify cs) of+ [x] -> return x+ _ -> error "totalsByLib: this is impossible" -classify :: [(Regex,String)] -> String -> String+-- --------------------------------------------------+-- | Parse clusters+readClusters :: FilePath -> IO [Cluster]+readClusters f = readFile f >>= return . rc . lines+ where rc [] = []+ rc [_] = error "odd number of cluster lines!"+ rc (c:ss:rest) = case head $ words c of+ ('>':name) -> (name,words ss) : rc rest+ _ -> error ("Cluster '"++c++"' didn't start with '>'")++classify :: LibTable -> String -> String classify ps str = case concatMap (class1 str) ps of [] -> error ("no match for "++str++" in library table")- [x] -> x++":"++str+ [x] -> x s@(_:_) -> error ("multiple matches for "++str++": "++show s) where class1 st (r,s) = maybe [] (const [s]) (matchRegex r st) -getClusters :: (String -> String) -> [String] -> [(String,[String])]-getClusters _ [] = []-getClusters _ [_] = error "odd number of cluster lines!"-getClusters f (c:ss:rest) = case head $ words c of- ('>':name) -> (name,map f $ words ss) : getClusters f rest- _ -> error ("Cluster '"++c++"' didn't start with '>'") --- transpose and print-writeClusters :: [Cluster] -> IO ()-writeClusters cs = do- let cs' = map (\(h,s) -> (h:s)) cs- cs'' = sortBy (\x y -> compare (length y) (length x)) cs'- unwords' = concat . intersperse ","- putStrLn $ unlines $ map unwords' $ transpose cs''+-- --------------------------------------------------+-- will need to match against all, to check for multiple matches+-- tag names with library+classClusters :: LibTable -> FilePath -> IO [Cluster]+classClusters ps f = return . map class1 =<< readClusters f+ where class1 (l,ss) = (l, map (\s -> classify ps s++":"++s) ss) -transpose :: [[String]] -> [[String]]-transpose [] = []-transpose ls' = let ls = filter (not . null) ls'- in map head ls : transpose (map tail ls)+-- | Output the clusters to stdout+writeTable :: [(String,[String])] -> IO ()+writeTable = putStrLn . unlines . map show1 + where+ show1 (name,stuff) = concat $ intersperse "\t" (name:stuff)++-- count by classification+countClusters :: LibTable -> [Cluster] -> [(String,[Int])]+countClusters lt = map count1+ where count1 (c,ss) = (c,map (\k -> findWithDefault 0 k (counts ss)) (map snd lt))+ counts ss = fromListWith (+) $ zip (map (classify lt) ss) $ repeat 1