diff --git a/README b/README --- a/README +++ b/README @@ -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 diff --git a/clustertools.cabal b/clustertools.cabal --- a/clustertools.cabal +++ b/clustertools.cabal @@ -1,9 +1,9 @@ Name: clustertools -Version: 0.1 +Version: 0.1.1 License: GPL License-file: LICENSE Author: Ketil Malde -Maintainer: Ketil Malde +Maintainer: Ketil Malde 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: . + * 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 @@ -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 diff --git a/src/Ace2Clusters.hs b/src/Ace2Clusters.hs new file mode 100644 --- /dev/null +++ b/src/Ace2Clusters.hs @@ -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 diff --git a/src/Ace2Fasta.hs b/src/Ace2Fasta.hs --- a/src/Ace2Fasta.hs +++ b/src/Ace2Fasta.hs @@ -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 diff --git a/src/AddSingletons.hs b/src/AddSingletons.hs --- a/src/AddSingletons.hs +++ b/src/AddSingletons.hs @@ -9,6 +9,7 @@ import Control.Monad import Data.Set +main :: IO () main = do (ref:fs) <- getArgs r <- readFile ref diff --git a/src/Cluscomp.lhs b/src/Cluscomp.lhs --- a/src/Cluscomp.lhs +++ b/src/Cluscomp.lhs @@ -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} diff --git a/src/ClusterLibs.hs b/src/ClusterLibs.hs --- a/src/ClusterLibs.hs +++ b/src/ClusterLibs.hs @@ -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 '>'") diff --git a/src/Filter.hs b/src/Filter.hs --- a/src/Filter.hs +++ b/src/Filter.hs @@ -1,5 +1,5 @@ -import System.SimpleArgs +import System.SimpleArgs (getArgs) import Data.Set hiding (filter,map,null) main :: IO () diff --git a/src/Xcerpt.lhs b/src/Xcerpt.lhs new file mode 100644 --- /dev/null +++ b/src/Xcerpt.lhs @@ -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 \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" + +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}