xml2x 0.4 → 0.4.1

raw patch · 2 files changed

+124/−3 lines, 2 files

Files

• src/Annotate.hs +121/−0
• xml2x.cabal +3/−3

@@ -0,0 +1,121 @@+{- | Link sequences (read from the xml file) to desired annotations +     (Currently support GO, todo support KEGG)+ -}+module Annotate where++import Data.ByteString.Lazy.Char8 (ByteString,unpack,pack)+import qualified Data.ByteString.Lazy.Char8 as B+import Data.Map  (Map)+import qualified Data.Map as M+import Data.Set  (Set)+import qualified Data.Set as S+import Data.List (nub)+import Data.Maybe++import Control.Monad+import Control.Arrow (second)++import Bio.Sequence.GeneOntology+import Bio.Sequence.KEGG+import qualified Bio.Alignment.BlastData as BD+import Bio.Alignment.BlastData (BlastResult)+import Bio.Alignment.BlastFlat +import Bio.Alignment.BlastXML as X -- new BlastFlat exports this too+import Data.Array.Unboxed+import Data.Word++import Options++type GoDefinitions = Map GoTerm (GoDef,[GoTerm])+type GoAnnotations = Map UniProtAcc (UArray Int Word16)+type KoAnnotations = Map UniProtAcc KO++{- | Read GO annotations using GOA.+   +   1. build a map of GO terms to description (25K lines - tiny)+   2. collect all proteins names from the XML file (1 pass over .xml)+   3. build a map of proteins to GO terms (1 pass over gene_assoc..), and add GO term desc+-}+read_go :: Options -> IO (GoDefinitions, GoAnnotations)+read_go opts  = do  -- (Opts (Just go) (Just goa),[xml],[])+  -- goterms are not necessary when we don't also use goanns+  case goann opts of+    Nothing -> return (M.empty,M.empty)+    Just goa -> do+      ts <- case godef opts of+              Nothing -> return M.empty+              Just obo -> do+                    t <- return . goTerms =<< (count opts) 1000 "reading GO terms: "+                           =<< readObo obo+                    when (isNothing $ M.lookup (GO 0) t) $ return () -- strictness hack+                    return t++      proteins <- getProts opts+      when (S.member (B.empty) proteins == False) $ return () -- strictness+      mt <- readGOA goa >>= (count opts) 1000 "reading GO annotations: "+            >>= return .  protTerms proteins+      when (isNothing $ M.lookup B.empty mt) $ return () -- strictness, again+      return (ts,mt)+  ++-- ----------------------------------------------------------------------+-- Build tables of annotation information+-- ----------------------------------------------------------------------++-- | Parse GO.terms_and_ids into a map+goTerms :: [(GoDef,[GoTerm])] -> GoDefinitions+goTerms = M.fromList . map go2keyval+    where go2keyval gd@(GoDef gt _ _,_) = (gt,gd)++-- | Read blast xml output and collect UniProt hits+getProts :: Options -> IO (Set UniProtAcc)+getProts opts = do+  hs <- mapM (\f -> (count opts) 100 ("reading accessions from '"++f++"': ")+                    . concatMap BD.hits . concatMap BD.results =<< X.readXML f) (inputs opts)+  return . S.fromList . map (chop . BD.subject) $ concat hs++-- | Convert "UniProt_ABCDEF Blah blah blah" to "ABCDEF"+--   Warning: only tested on this format, may break!+chop :: ByteString -> ByteString+chop = B.copy . B.drop 1 . B.dropWhile (/='_') . head . B.words++-- | Using GO definitions and the set of uniprot terms, scan the associations file and+--   snarf all relevant data.  Requires the annotations to be grouped by protein.+protTerms :: (Set UniProtAcc) -> [Annotation] -> GoAnnotations+protTerms ps as = M.fromAscList . map toArray . partitions . filter isMember $ as+    where toArray (Ann up (GO gt1) _:xs) = let gs = nub (gt1:[ gt | Ann _ (GO gt) _ <- xs ])+                                               a  = listArray (0,length gs-1) $ map fromIntegral gs+                                           in a `seq` (pack . unpack $ up,a)+          toArray [] = error "cannot read an empty list of annotations!"+          isMember (Ann up _ _) = up `S.member` ps+          partitions xs@(Ann x _ _:_) = let (p1,pps) = span (\(Ann y _ _) -> y==x) xs+                                        in p1:partitions pps+          partitions [] = []+++{- | Read KEGG annotations +   KEGG is specified as abbreviated organism prefixes. +   Files to be read are $prefix_ko.list and $prefix_uniprot.list.+ -}+read_kegg :: Options -> IO KoAnnotations+read_kegg opts = M.unions `fmap` mapM build_kegg (ko opts)+    where -- build KEGG info for a single organism+          build_kegg org = do +            k <- return . M.fromList =<<+                 (count opts) 1000 ("parsing KEGG information for '"++org++"'.. ")+                     =<< merge_files org+            when (isNothing $ M.lookup (B.pack "") k) $ return () -- strictness hack+            return k+          -- parse and merge files; they seem to be lexicographically sorted+          merge_files org = do+            g2u <- map (second decodeUP) `fmap` genReadKegg (org++"_uniprot.list") +            g2k <- map (second decodeKO) `fmap` genReadKegg (org++"_ko.list")+            return $ merge g2u g2k++          -- merge sort on fst, returning snd            +          merge :: [(ByteString,UniProtAcc)] -> [(ByteString,KO)] -> [(UniProtAcc,KO)]+          merge (x:xs) (y:ys) | fst x < fst y = merge xs (y:ys)+                              | fst x > fst y = merge (x:xs) ys+                              | otherwise     = (snd x,snd y) : merge xs ys+          merge [] _ = []+          merge _ [] = []

@@ -1,5 +1,5 @@ Name:           xml2x-Version:        0.4+Version:        0.4.1 License:        GPL License-File:   LICENSE @@ -26,5 +26,5 @@ Main-Is:        Xml2X.hs Extensions:     Rank2Types Hs-Source-Dirs: src-Other-modules:  Html, Options, Tabulate-Ghc-Options:    -Wall -O2 -smp+Other-modules:  Html, Options, Tabulate, Annotate+Ghc-Options:    -Wall -O2 -threaded