clustertools 0.1.2 → 0.1.5
raw patch · 6 files changed
+151/−94 lines, 6 filesdep −haskell98dep ~basedep ~bio
Dependencies removed: haskell98
Dependency ranges changed: base, bio
Files
- clustertools.cabal +18/−12
- src/Ace2Fasta.hs +2/−3
- src/ClusterLibs.hs +7/−58
- src/Formats.hs +57/−0
- src/Statistics.hs +39/−0
- src/Xcerpt.lhs +28/−21
clustertools.cabal view
@@ -1,5 +1,5 @@ Name: clustertools-Version: 0.1.2+Version: 0.1.5 License: GPL License-file: LICENSE Author: Ketil Malde@@ -22,18 +22,17 @@ . * ace2fasta - parse an ACE assembly, and output each assembly in a separate FASTA file .- * ace2clusters - parse an ACE assembly, and output clusters in TGICL format- .+ * ace2clusters - parse an ACE assembly, and output clusters in TGICL format+ . * clusterlibs - given a table of regular expressions and library names, along with a- clustering (TGICL-format), output a table of clusters with the library name- prepended to the sequences.+ clustering (TGICL-format), output a table of cluster sizes per library. .- * xcerpt - extract sequences from a list of sequence labels.- .+ * xcerpt - extract sequences from a list of sequence labels.+ . The Darcs repository is at: <http://malde.org/~ketil/biohaskell/cluster_tools>. Homepage: http://malde.org/~ketil/ -Tested-With: GHC==6.8.2+Tested-With: GHC==6.10.4 Build-Type: Simple Cabal-Version: >= 1.2 @@ -42,39 +41,46 @@ Executable filter Hs-source-dirs: src Main-Is: Filter.hs- Build-Depends: base>3, containers, simpleargs>=0.1- ghc-options: -Wall -O2 -funbox-strict-fields+ Build-Depends: base >= 4 && < 5, containers, simpleargs>=0.1+ ghc-options: -Wall -funbox-strict-fields Executable clusc Hs-source-dirs: src Main-Is: Cluscomp.lhs Build-Depends: bytestring+ Ghc-Options: -Wall Executable add_single Hs-source-dirs: src Main-Is: AddSingletons.hs+ Ghc-Options: -Wall Executable ace2contigs Hs-source-dirs: src Main-Is: Ace2Contigs.hs Build-depends: bio>=0.3.3.4+ Ghc-Options: -Wall Executable ace2fasta Hs-source-dirs: src Main-Is: Ace2Fasta.hs- Build-depends: bio>=0.3.3.4+ Build-depends: bio>=0.4+ Ghc-Options: -Wall Executable ace2clusters Hs-source-dirs: src Main-Is: Ace2Clusters.hs Build-depends: bio>=0.3.3.4+ Ghc-Options: -Wall Executable clusterlibs Hs-source-dirs: src Main-Is: ClusterLibs.hs+ Other-modules: Statistics, Formats Build-depends: regex-compat, QuickCheck+ Ghc-Options: -Wall Executable xcerpt Hs-source-dirs: src Main-Is: Xcerpt.lhs- Build-depends: haskell98+ Ghc-Options: -Wall
src/Ace2Fasta.hs view
@@ -9,19 +9,18 @@ import Bio.Sequence import qualified Data.ByteString.Lazy.Char8 as B import System.SimpleArgs (getArgs)-import System.IO main :: IO () main = do a <- readACE =<< getArgs mapM_ (\ss -> writeFasta (B.unpack $ seqlabel $ head ss) ss) (map pad_all $ concat a) -pad_all :: Assembly -> [Sequence]+pad_all :: Assembly -> [Sequence Nuc] pad_all asm = pad 21 (contig asm) : map padRead (fragments asm) where padRead (off,_dir,sq,gs) = pad (20+off) (sq,gs) -- note: they are already complemented -- generate a sequence with '-' for gaps-pad :: Offset -> (Sequence,Gaps) -> Sequence+pad :: Offset -> (Sequence Nuc,Gaps) -> Sequence Nuc pad off (Seq n s _,gs) = let s' = (if off < 0 then B.drop (negate off) else B.append (B.replicate off '-')) (insertGaps '-' (s,gs)) in (Seq n s' Nothing)
src/ClusterLibs.hs view
@@ -2,16 +2,15 @@ -- calculate clusters by library, using the lib table -- read patterns from a table, first line is header +module Main where+ import System.Environment (getArgs)-import Text.Regex-import Data.List (elemIndex,sortBy,intersperse)-import Data.Map hiding (map)+import Data.List (intersperse) import Numeric (showFFloat) import Statistics--type LibTable = [(Regex,String)]-type Cluster = (String,[String])+import Formats (LibTable, Cluster+ , readPatternTable, totalsByLib, readClusters, countClusters, classify) main :: IO () main = do@@ -23,7 +22,8 @@ writeTable [("","sum":map snd pat++["significance"])] writeTable . map (decorate counts) . countClusters pat =<< readClusters cs --- counts is +-- | Using count by library, takes a cluster, calculates significance scores +-- and formats it for output. decorate :: (String,[Int]) -> (String,[Int]) -> (String,[String]) decorate counts (x,ys) = let fractions = [ fromIntegral x / fromIntegral (sum $ snd counts) | x <- snd counts]@@ -32,11 +32,6 @@ isSignificant obs frac = showFFloat (Just 3) (pvalue (1-frac) clustersize (clustersize-obs)) "" in (x,map show (clustersize:ys) ++ [significance]) -{-- clus <- classClusters pat cs- writeTable clus--}- -- -------------------------------------------------- -- | Calculate the p value of observing k sequences -- from a library, given an a priori background distribution. @@ -44,47 +39,6 @@ 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- (Just x,Just y) -> (x,y)- _ -> error ("Need both 'Pattern' and 'Name' headers in '"++f++"'")- z l | length l < max p n = error ("Line in library table too short:\n"++show l)- | otherwise = (mkRegex (l!!p),l!!n)- return $ map z ls---- ----------------------------------------------------- | 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"---- ----------------------------------------------------- | 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- s@(_:_) -> error ("multiple matches for "++str++": "++show s)- where class1 st (r,s) = maybe [] (const [s]) (matchRegex r st)----- -------------------------------------------------- -- will need to match against all, to check for multiple matches -- tag names with library classClusters :: LibTable -> FilePath -> IO [Cluster]@@ -97,8 +51,3 @@ 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
+ src/Formats.hs view
@@ -0,0 +1,57 @@+-- | Support for reading various data +module Formats where+ +import Text.Regex+import Data.List (elemIndex)+import Data.Map hiding (map)++type LibTable = [(Regex,String)] -- | pattern and name+type Cluster = (String,[String]) -- | name and list of reads++-- --------------------------------------------------+-- | Classify a read according to the LibTable+classify :: LibTable -> String -> String+classify ps str = case concatMap (class1 str) ps of+ [] -> error ("no match for "++str++" in library table")+ [x] -> x+ s@(_:_) -> error ("multiple matches for "++str++": "++show s)+ where class1 st (r,s) = maybe [] (const [s]) (matchRegex r st)++-- --------------------------------------------------+-- | 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+ (Just x,Just y) -> (x,y)+ _ -> error ("Need both 'Pattern' and 'Name' headers in '"++f++"'")+ z l | length l < max p n = error ("Line in library table too short:\n"++show l)+ | otherwise = (mkRegex (l!!p),l!!n)+ return $ map z ls++-- --------------------------------------------------+-- | 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"++-- 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++-- --------------------------------------------------+-- | 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 '>'")
+ src/Statistics.hs view
@@ -0,0 +1,39 @@+-- | Statistics for calculating probabilities from the binomial distribution++module Statistics where+import Test.QuickCheck hiding (choose)++-- | Number of combinations for choosing k from a set of n objects. Needs to use Integer to+-- avoid running out of Int bits.+choose :: Integral i => i -> i -> Integer+n' `choose` k' = let (n,k) = (fromIntegral n', fromIntegral k')+ in product [k+1..n] `div` product [1..n-k]++-- | Calculate binomial probability of observing exactly k positives from n trials,+-- with a probability p for a positive result in each trial.+binomial :: Integral i => Double -> i -> i -> Double+binomial p n k -- | p*fromIntegral n > 20 && (1-p)*fromIntegral n > 20 = normal_apporx p n k+ | k <= n = fromIntegral (n `choose` k) * (p**fromIntegral k) * ((1-p)** fromIntegral (n-k))+ | otherwise = error ("binomial: can't observe more than 'n' positives:"++show n++" "++show k)++prop_binom p n k = k <= n && p>=0 && p <= 1 ==> sum [ binomial p n k | k <- [0..n]] - 1 < 0.001+++-- | Calculate cumulative binomial probability of achieving k or fewer positive results.+cumbin :: Integral i => Double -> i -> i -> Double+cumbin p n k = sum [binomial p n x | x <- [0..k]]++prop_cumbin p n k = k <= n && n-k-1 <= n && p>=0 && p <= 1+ ==> cumbin p n k + cumbin (1-p) n (n-k-1) - 1 < 0.001++{- crap! this won't work, you can't just sum normals, goddamit.+normal_approx p n k = pdf_normal (fromIntegral n*p) (fromIntegral n*p*(1-p)) (fromIntegral k)+pdf_normal mu sigma x = exp(negate((x-mu)^2/(2*sigma^2)))/(sigma*sqrt(2*pi))++invcumnorm mu sigma z = mu + search (-limit*sigma) (limit*sigma)+ where search a b = let c = (a+b)/2+ cn = cumnorm 0 sigma c+ in if abs (z - cn) < 10*epsilon || abs (a-b) < epsilon then c+ else if cn > z then search a c+ else search c b+-}
src/Xcerpt.lhs view
@@ -5,48 +5,55 @@ module Main where -import Data.Set hiding (null,filter)-import System-import Data.List (foldl')-import IO+import Data.Set hiding (null,filter,partition)+import System.Environment (getArgs)+import Data.List (foldl',partition)+import System.IO usage :: String-usage = "xcerpt <labels> <input>\n\n" +++usage = "xcerpt [-r|-f] [<labels>] <input>\n\n" ++ "where <labels> is the name of a file containing the IDs\n" ++ "of sequences to extract from the <input> file.\n" ++- "The result ends up in <input>.match and <input>.rest"+ "The default is to output matching sequences, -r outputs the non-matching sequences" +++ "With -f, the result ends up in xcerpt.match and xcerpt.rest instead of standard output." main :: IO () main = do- args <- System.getArgs- if length args /= 2 then - putStrLn usage- else do- d <- readFile (args!!0)- let dict = mkdict d- xcerpt dict (args!!1)+ (opts,args) <- return . partition (\a->not (null a) && head a=='-') =<< getArgs+ let dict = case length args of + 1 -> hGetContents stdin+ 2 -> readFile (args!!0)+ _ -> error usage+ d <- return . mkdict =<< dict+ (match_out,rest_out) <- case opts of + [] -> return (Just stdout,Nothing)+ ["-r"] -> return (Nothing,Just stdout)+ ["-f"] -> do m <- openFile "xcerpt.match" WriteMode + r <- openFile "xcerpt.rest" WriteMode+ return (Just m, Just r)+ _ -> error usage+ xcerpt d match_out rest_out (last args)+ maybe (return ()) hClose match_out+ maybe (return ()) hClose rest_out mkdict :: String -> Set String mkdict = foldl' (flip insert) empty . words -xcerpt :: Set String -> String -> IO ()-xcerpt dict input = do- m <- openFile (input++".match") WriteMode- r <- openFile (input++".rest") WriteMode+xcerpt :: Set String -> Maybe Handle -> Maybe Handle -> String -> IO ()+xcerpt dict m r input = do i <- readFile input xtr dict m r $ filter (\l->(not.null) l && head l /= '#') $ lines i -xtr _ m r [] = do- hClose m- hClose r+xtr _ m r [] = return () xtr _ _ _ [x] = error ("Odd number of lines?\n"++x) xtr d m r (l1:ls) = if head l1 == '>' then let f = if (drop 1 $ head $ words l1) `member` d then m else r in do let (sequence,rest) = break ((=='>').head) ls- hPutStr f $ unlines (l1:sequence)+ case f of Just x -> hPutStr x $ unlines (l1:sequence)+ Nothing -> return () xtr d m r rest else error ("Not a FASTA header:\n"++l1)