diff --git a/clustertools.cabal b/clustertools.cabal --- a/clustertools.cabal +++ b/clustertools.cabal @@ -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: . 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 diff --git a/src/Ace2Fasta.hs b/src/Ace2Fasta.hs --- a/src/Ace2Fasta.hs +++ b/src/Ace2Fasta.hs @@ -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) diff --git a/src/ClusterLibs.hs b/src/ClusterLibs.hs --- a/src/ClusterLibs.hs +++ b/src/ClusterLibs.hs @@ -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 diff --git a/src/Formats.hs b/src/Formats.hs new file mode 100644 --- /dev/null +++ b/src/Formats.hs @@ -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 '>'") diff --git a/src/Statistics.hs b/src/Statistics.hs new file mode 100644 --- /dev/null +++ b/src/Statistics.hs @@ -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 +-} diff --git a/src/Xcerpt.lhs b/src/Xcerpt.lhs --- a/src/Xcerpt.lhs +++ b/src/Xcerpt.lhs @@ -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 \n\n" ++ +usage = "xcerpt [-r|-f] [] \n\n" ++ "where is the name of a file containing the IDs\n" ++ "of sequences to extract from the file.\n" ++ - "The result ends up in .match and .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)