clustertools 0.1 → 0.1.1
raw patch · 9 files changed
+124/−21 lines, 9 filesdep +haskell98dep ~bionew-component:exe:ace2clustersnew-component:exe:xcerptnew-uploader
Dependencies added: haskell98
Dependency ranges changed: bio
Files
- README +7/−1
- clustertools.cabal +19/−4
- src/Ace2Clusters.hs +17/−0
- src/Ace2Fasta.hs +2/−2
- src/AddSingletons.hs +1/−0
- src/Cluscomp.lhs +12/−7
- src/ClusterLibs.hs +12/−6
- src/Filter.hs +1/−1
- src/Xcerpt.lhs +53/−0
README view
@@ -31,8 +31,14 @@ and the corresponding quality information (.qual). ace2fasta - parse an ACE assembly, and output each assembly in a separate- fasta formatted file, with the necessary gaps inserted to align the+ FASTA formatted file, with the necessary gaps inserted to align the sequences (suitable for import into e.g. Seaview)++ace2clusters - parse an ACE assembly, and output clusters composed of the+ sequences used for each contig. The format is similar to TGICL's, + with cluster output as one line consisting of a '>' and the contig name,+ and the next line containing the names of the sequences that comprise+ the cluster. clusterlibs - given a table of regular expressions and library names, along with a clustering (TGICL-format), output a table of clusters
clustertools.cabal view
@@ -1,9 +1,9 @@ Name: clustertools-Version: 0.1+Version: 0.1.1 License: GPL License-file: LICENSE Author: Ketil Malde-Maintainer: Ketil Malde <ketil@ii.uib.no>+Maintainer: Ketil Malde <ketil@malde.org> Category: Bioinformatics Synopsis: Tools for manipulating sequence clusters@@ -22,11 +22,15 @@ . * 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+ . * 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. .- The Darcs repository is at: <http://malde.org/~ketil/cluster_tools>.+ * 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@@ -53,13 +57,24 @@ Executable ace2contigs Hs-source-dirs: src Main-Is: Ace2Contigs.hs- Build-depends: bio>=0.3.3+ Build-depends: bio>=0.3.3.4 Executable ace2fasta Hs-source-dirs: src Main-Is: Ace2Fasta.hs+ Build-depends: bio>=0.3.3.4 +Executable ace2clusters+ Hs-source-dirs: src+ Main-Is: Ace2Clusters.hs+ Build-depends: bio>=0.3.3.4+ Executable clusterlibs Hs-source-dirs: src Main-Is: ClusterLibs.hs Build-depends: regex-compat++Executable xcerpt+ Hs-source-dirs: src+ Main-Is: Xcerpt.lhs+ Build-depends: haskell98
+ src/Ace2Clusters.hs view
@@ -0,0 +1,17 @@+-- parse ACE files, and generate clusters á la TGICL, but based on contigs++module Main where++import Bio.Alignment.ACE+import Bio.Sequence+import System.SimpleArgs (getArgs)++main :: IO ()+main = getArgs >>= readACE >>= putStr . concatMap (concatMap decluster)++decluster :: Assembly -> String+decluster a = ">"++(label $ fst $ contig a)++"\n"+ ++ unwords (map (label . thrd) $ fragments a)++"\n"+ where+ thrd (_,_,s,_) = s+ label = toStr . seqlabel
src/Ace2Fasta.hs view
@@ -3,8 +3,8 @@ -- (for importing into seaview) module Main where-import Prelude hiding (reads) +import Bio.Alignment.AlignData (Gaps,insertGaps) import Bio.Alignment.ACE import Bio.Sequence import qualified Data.ByteString.Lazy.Char8 as B@@ -17,7 +17,7 @@ mapM_ (\ss -> writeFasta (B.unpack $ seqlabel $ head ss) ss) (map pad_all $ concat a) pad_all :: Assembly -> [Sequence]-pad_all asm = pad 21 (contig asm) : map padRead (reads asm)+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
src/AddSingletons.hs view
@@ -9,6 +9,7 @@ import Control.Monad import Data.Set +main :: IO () main = do (ref:fs) <- getArgs r <- readFile ref
src/Cluscomp.lhs view
@@ -31,6 +31,7 @@ putStrLn (pad n "clustering" ++ "\tVI\tMIF\tJaccard\tRand\tFowMal\tSkew") mapM_ (comp n ref) args +pad :: Int -> String -> String pad i s = s ++ replicate (i-length s) ' ' comp :: Int -> Clustering -> FilePath -> IO ()@@ -98,14 +99,14 @@ type Pairs = (Double,Double,Double,Double) -- a b c, and d -jaccard, rand, fowmal, skew :: Pairs -> Double+jaccard, rand, fowmal, hubert, skew :: Pairs -> Double jaccard (a,b,c,_) = a/(a+b+c) rand (a,b,c,d) = (a+d)/(a+b+c+d)-fowmal (a,b,c,d) = a/sqrt((a+b)*(a+c)) -- Fowlkes-Mallows index+fowmal (a,b,c,_) = a/sqrt((a+b)*(a+c)) -- Fowlkes-Mallows index -- this is the same as fowmal (why?)! hubert (a,b,c,d) = (a*d - b*c)/sqrt((a+b)*(c+d)*(a+c)*(b+d)) -- Gamma stat? -skew (a,b,c,d) = (b-c)/sqrt(b+c)/sqrt d+skew (_,b,c,d) = (b-c)/sqrt(b+c)/sqrt d mkPairs :: Clustering -> Matrix -> Pairs mkPairs k m = let (a,b,c,ct) = foldl (jc1 k) (0,0,0,0) m@@ -116,6 +117,8 @@ jc1 k j row = let row' = filter (isJust.fst) row l = sum $ map (len.snd) row' in foldr (jc2 k l) j row'++jc2 :: Clustering -> Double -> (Maybe CID,[Sequence]) -> Pairs -> Pairs jc2 k lc (Just i,ss) (a,b,c,ct) = (a+ls*(ls-1),b+(ls*(lk-ls)),c+(ls*(lc-ls)),ct+ls) where ls = len ss lk = len $ fromJust (seqs k i) -- fixme! only seqs found in C!@@ -155,20 +158,22 @@ let r' = map (len.snd) $ filter (isJust.fst) r in iscores (accum m r) (n+sum r',hc,hk+mlog (sum r'),hck+sum (map mlog r')) rs-iscores m (n,hc,hk,hck) [] =+iscores m (n,_hc,hk,hck) [] = (n,(sum $ map mlog $ elems m),hk,hck) -- accumulate the column sums for the matrix. This must be done -- explicitly, since one clustering may contain sequences not present -- in the other accum :: Map Int Double -> Row -> Map Int Double-accum m row = foldr accum1 m row- where accum1 (Just cid,seqs) m = insertWith (+) cid (len seqs) m- accum1 (Nothing,_) m = m+accum m' row = foldr accum1 m' row+ where accum1 (Just c,ss) m = insertWith (+) c (len ss) m+ accum1 (Nothing,_) m = m +mlog :: Double -> Double mlog 0 = 0 -- error "log 0!" mlog x = x*log x/log 2 +len :: [a] -> Double len = fromIntegral . length \end{code}
src/ClusterLibs.hs view
@@ -20,8 +20,11 @@ type Cluster = (String,[String]) +main :: IO () main = do- [ps,cs] <- getArgs `catch` error "Usage: clusterlibs ps cs"+ [ps,cs] <- do [xs,ys] <- getArgs + return [xs,ys] + `catch` error "Usage: clusterlibs ps cs" pat <- readPatternTable ps clus <- classClusters pat cs writeClusters clus@@ -32,7 +35,8 @@ 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 = (mkRegex (l!!p),l!!n)+ 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 -- will need to match against all, to check for multiple matches@@ -40,14 +44,16 @@ classClusters :: [(Regex,String)] -> FilePath -> IO [Cluster] classClusters ps f = return . getClusters (classify ps) . lines =<< readFile f +classify :: [(Regex,String)] -> String -> String classify ps str = case concatMap (class1 str) ps of- [] -> error ("no match for "++str)+ [] -> error ("no match for "++str++" in library table") [x] -> x++":"++str s@(_:_) -> error ("multiple matches for "++str++": "++show s)- where class1 str (r,s) = maybe [] (const [s]) (matchRegex r str)+ where class1 st (r,s) = maybe [] (const [s]) (matchRegex r st) -getClusters f [] = []-getClusters f [_] = error "odd number of cluster lines!"+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 '>'")
src/Filter.hs view
@@ -1,5 +1,5 @@ -import System.SimpleArgs+import System.SimpleArgs (getArgs) import Data.Set hiding (filter,map,null) main :: IO ()
+ src/Xcerpt.lhs view
@@ -0,0 +1,53 @@++Extract sequences matching a list of labels++\begin{code}++module Main where++import Data.Set hiding (null,filter)+import System+import Data.List (foldl')+import IO++usage :: String+usage = "xcerpt <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"++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)++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+ 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 _ _ _ [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)+ xtr d m r rest+ else error ("Not a FASTA header:\n"++l1)++\end{code}