bio 0.3.3.2 → 0.3.3.4

raw patch · 15 files changed

+403/−133 lines, 15 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• Bio/Alignment/ACE.hs +27/−19
• Bio/Alignment/Blast.hs +6/−1
• Bio/Alignment/BlastFlat.hs +6/−0
• Bio/Alignment/BlastXML.hs +12/−2
• Bio/Alignment/Matrices.hs +21/−6
• Bio/Alignment/QAlign.hs +74/−10
• Bio/Sequence.hs +1/−0
• Bio/Sequence/Fasta.hs +2/−2
• Bio/Sequence/GOA.hs +4/−85
• Bio/Sequence/GeneOntology.hs +208/−0
• Bio/Sequence/KEGG.hs +18/−0
• Bio/Sequence/SeqData.hs +13/−0
• Bio/Util.hs +7/−6
• Bio/Util/Parsex.hs +1/−1
• bio.cabal +3/−1

@@ -9,19 +9,22 @@     Briefly: each field is a line starting with a two letter code,             in some cases followed by data lines termintated by a blank line.-          -- AS contigs reads-          --  CO contig_name bases reads segments compl (CAP3: segments=0)-          --    sequence-          --  BQ-          --    base_qualities-          --  AF read1 compl padded_start_consensus (negatives meaning?)-          --  AF read2 ..-          --  BS segments-          --  RD read1 bases info_items info_tags (latter two set to 0 by CAP3)-          --    sequence-          --  QA read1 qual_start qual_end align_start align_end-          --  DS (phred header? left empty by CAP3)-          --  RD read2 ...+   Here's an brief example how an ACE file looks like:++   @+          AS contigs reads+          CO contig_name bases reads segments compl (CAP3: segments=0)+          sequence+          BQ base_qualities+          AF read1 compl padded_start_consensus (negatives meaning?)+          AF read2 ..+          BS segments+          RD read1 bases info_items info_tags (latter two set to 0 by CAP3)+          sequence+          QA read1 qual_start qual_end align_start align_end+          DS (phred header? left empty by CAP3)+          RD read2 ...+   @ -}  {-# LANGUAGE CPP #-}@@ -124,11 +127,14 @@   sd <- sdata   let (sd',gaps) = extractGaps sd   many blank-  bq-  sq <- qdata+  -- Vector NTI produces ACE without BQ for the contig+  msq <- do bq+            sq <- qdata+            return (Just sq)+       <|> return Nothing   many blank   -- todo: gaps?-  return (Seq name sd' (Just sq),gaps)+  return (Seq name sd' msq,gaps)  co, sdata, qdata :: AceParser Str co = parse1 (\t -> case t of CO name a b c _compl -> do@@ -205,6 +211,8 @@   return (Seq rn s Nothing,gaps) -- huh?  -- | parse each read (RD, QA, DS)+--   Vector NTI appears to insert solitary RDs, sometimes even without any sequence data!?+--   This is not supported at this point. rd :: AceParser (Str,Int,Int,Int) rd = parse1 (\t -> case t of RD a b c d -> do [x,y,z] <- readInts [b,c,d]                                               return (a,x,y,z)@@ -224,9 +232,9 @@ -- Tokenise a single line -- todo: error on incorrect (partial) format, error reports with line number tokenize1 :: ByteString -> ACE-tokenize1 l | B.null l = Empty-            | otherwise = let (h:ws) = words l-                          in case (unpack h,ws) of+tokenize1 l = case words l of+                [] -> Empty+                (h:ws) -> case (unpack h,ws) of                             ("AS",[cs,rs])               -> AS cs rs                             ("CO",[nm,bss,rs,segs,comp]) -> CO nm bss rs segs comp                             ("BQ",[])                    -> BQ

@@ -36,6 +36,7 @@ import Bio.Util (splitWhen)  -- String constant used in parsing+str_query, str_gt, str_score, str_refer, str_datab, str_search :: ByteString str_query = B.pack "Query= " str_gt    = B.pack ">" str_score = B.pack " Score = "@@ -43,7 +44,7 @@ str_datab = B.pack "Database: " str_search= B.pack "Searching" -#define BUG(C_,M_) (error ("Program error - '"++C_++"' failed: "++M_++". Location: "++__FILE__++" line: "++ show __LINE__))+#define BUG(C_,M_) (error ("Program error - '"++C_++"' failed: "++M_++". Location: "++__FILE__++" line: "++ show (__LINE__::Int))) #define readX (\s -> case [ x | (x,t) <- reads (B.unpack s), ("","") <- lex t] of { [x] -> x ; [] -> BUG("read",("no parse for "++B.unpack s)); _ -> BUG("read","ambigous parse")}) -------------------------------------------------- -- Splitting up the input in its relevant parts --@@ -71,6 +72,7 @@   -- parse metadata from the preamble+parse_preamble :: [ByteString] -> BlastResult parse_preamble (l:ls) = BlastResult { blastprogram = w1, blastversion = w2, blastdate = w3                                     , blastreferences = B.concat $ map (B.drop (B.length str_refer)) rs                                     , database = B.drop (B.length str_datab) d@@ -81,17 +83,20 @@           d = case filter (isPrefixOf str_datab) records of [d1] -> d1; _ -> B.pack "<unknown>"  -- parse the blast record for one query sequence+parse_query :: [ByteString] -> BlastRecord parse_query ls = let ((s1:s2:_):hs) = qhits ls                  in BlastRecord { query = B.drop (B.length str_query) s1                                 , qlength = readX $ B.takeWhile isDigit $ (!!1) $ B.split '(' s2                                 , hits = map parse_hit hs }  -- parse a hit against a sequence in the data base+parse_hit :: [ByteString] -> BlastHit parse_hit ls = let (h:ms) = hmatches ls                in BlastHit { subject = B.concat h, slength = -1, matches = map parse_match ms}  -- parse a match between a query and a subject -- the format is a bit variable, it seems+parse_match :: [ByteString] -> BlastMatch parse_match (l1:l2:l3:_) = let     hdr = words $ B.concat [l1,l2,l3]     [bs,ev,ident,str1,str2] = map (hdr!!) [2,7,10] ++ map (reverse hdr!!) [2,0] -- Strand is last, ignore optional Gap

@@ -17,6 +17,8 @@     (      -- * The BlastFlat data type       BlastFlat(..)+    -- * Read XML format+    , readXML     -- * Convert from hierarchical to flat structure     , flatten     -- * Re-exports from the hierarchical module ("Bio.Alignment.BlastData")@@ -27,6 +29,7 @@     )where  import qualified Bio.Alignment.BlastData as B+import qualified Bio.Alignment.BlastXML as X import Data.ByteString.Lazy.Char8 (empty)  -- | The BlastFlat data structure contains information about a single match@@ -38,6 +41,9 @@     , q_from, q_to, h_from, h_to :: !Int     , aux :: !B.Aux         }++readXML :: FilePath -> IO [BlastFlat]+readXML f = return . concatMap (flatten . B.results) =<< X.readXML f   -- | Convert BlastRecords into BlastFlats (representing a depth-first traversal of the  --   BlastRecord structure.)

@@ -1,5 +1,15 @@--- Parse blast XML output using tagsoup+{- | Parse blast XML output. +   If you use a recent version of NCBI BLAST and specify XML output (blastall -m 7),+   this module should be able to parse the result into a hierarchical 'BlastResult'+   structure.++   While the process may consume a bit of memory, the parsing is lazy,+   and file sizes of several gigabytes can be parsed (see e.g. the+   xml2x tool for an example).  To parse XML, we use+   'Text.HTML.TagSoup'. +-}+ module Bio.Alignment.BlastXML (readXML) where  import Bio.Alignment.BlastData@@ -7,9 +17,9 @@ import qualified Data.ByteString.Lazy.Char8 as B import Text.HTML.TagSoup import Control.Monad- import Control.Parallel +-- | Parse BLAST results in XML format readXML :: FilePath -> IO [BlastResult] readXML fp = do      ts <- return . parseTags =<< readFile fp

@@ -4,28 +4,39 @@    When in doubt, use BLOSUM62.    Consult <http://www.ncbi.nlm.nih.gov/blast/blast_whatsnew.shtml#20051206>    for some hints on good parameters for nucleotide alignments.+ +   See also <http://en.wikipedia.org/wiki/Substitution_matrix> for a+   summary about the difference between the different matrices.  -}  module Bio.Alignment.Matrices     (-    -- * BLOSUM matrices (Henikoff and Henikoff)+    -- * BLOSUM matrices +    -- | For BLOSUM matrices, the associated number determines the+    -- similarity of the sequences the matrices are derived from.+    --+    -- Henikoff, S. and Henikoff, J. Amino acid substitution matrices from protein blocks. Proc. Natl. Acad. Sci. USA. 89(biochemistry): 10915 - 10919 (1992).     blosum45, blosum62, blosum80 -    -- * PAM matrices (Dayhoff et al)+    -- * PAM matrices +    -- | For PAM matrics, the number indicates the number of mutations+    --   that have occurred between the sequences that are compared.+    --+    -- Dayhoff, M.O., Schwartz, R.M., Orcutt, B.C. A model of evolutionary change in proteins. In \"Atlas of Protein Sequence and Structure\" 5(3) M.O. Dayhoff (ed.), 345 - 352 (1978).     , pam30, pam70      -- * BLASTn defaults, for nucleotide sequences     , blastn_default -    -- * Generic and simple "matrix" generator+    -- * Generic and simple matrix generator     , simpleMx     ) where  import Bio.Alignment.AlignData (Chr) import qualified Data.Map as M --- | The standard BLOSUM45 matrix.-blosum45, blosum62, blosum80, pam30, pam70 :: (Char,Char) -> Int+-- | BLOSUM45 matrix, suitable for distantly related sequences+blosum45 :: (Char,Char) -> Int blosum45 m = M.findWithDefault (-5) m $ M.fromList [(('A','A'),5),     (('A','B'),-1), (('A','C'),-1), (('A','D'),-2), (('A','E'),-1),     (('A','F'),-2), (('A','G'),0), (('A','H'),-2), (('A','I'),-1),@@ -161,6 +172,7 @@     (('Z','W'),-2), (('Z','X'),-1), (('Z','Y'),-2), (('Z','Z'),4)]  -- | The standard BLOSUM62 matrix.+blosum62 :: (Char,Char) -> Int blosum62 m = M.findWithDefault (-4) m $ M.fromList      [(('A','A'),4), (('A','B'),-2), (('A','C'),0), (('A','D'),-2),      (('A','E'),-1), (('A','F'),-2), (('A','G'),0), (('A','H'),-2),@@ -296,7 +308,8 @@      (('Z','V'),-2), (('Z','W'),-3), (('Z','X'),-1), (('Z','Y'),-2),      (('Z','Z'),4)] --- | The standard BLOSUM80 matrix.+-- | BLOSUM80 matrix, suitable for closely related sequences. +blosum80 :: (Char,Char) -> Int blosum80 m = M.findWithDefault (-6) m $ M.fromList [(('A','A'),5),     (('A','B'),-2), (('A','C'),-1), (('A','D'),-2), (('A','E'),-1),     (('A','F'),-3), (('A','G'),0), (('A','H'),-2), (('A','I'),-2),@@ -432,6 +445,7 @@     (('Z','W'),-4), (('Z','X'),-1), (('Z','Y'),-3), (('Z','Z'),4)]  -- | The standard PAM30 matrix+pam30 :: (Char,Char) -> Int pam30 m = M.findWithDefault (-17) m $ M.fromList [(('A','A'),6),     (('A','B'),-3), (('A','C'),-6), (('A','D'),-3), (('A','E'),-2),     (('A','F'),-8), (('A','G'),-2), (('A','H'),-7), (('A','I'),-5),@@ -568,6 +582,7 @@     (('Z','Z'),6)]  -- | The standard PAM70 matrix.+pam70 :: (Char,Char) -> Int pam70 m = M.findWithDefault (-11) m $ M.fromList [(('A','A'),5),     (('A','B'),-1), (('A','C'),-4), (('A','D'),-1), (('A','E'),-1),     (('A','F'),-6), (('A','G'),0), (('A','H'),-4), (('A','I'),-2),

@@ -8,6 +8,8 @@     This module performs sequences alignments, takes quality values into    account.++   See also <http://bioinformatics.oxfordjournals.org/cgi/content/abstract/btn052v1>. -} {-# LANGUAGE CPP, ParallelListComp #-} {-# OPTIONS_GHC -fexcess-precision #-}@@ -15,18 +17,28 @@ -- #define DEBUG  module Bio.Alignment.QAlign (-   -- * Smith-Waterman, or locally optimal alignment with affine gaps+   -- * Smith-Waterman+   -- | Locally optimal alignment with affine gaps, i.e. best infix match.    local_score, local_align -   -- * Needleman-Wunsch, or globally optimal alignment with affine gaps+   -- * Needleman-Wunsch +   -- | Globally optimal alignment with affine gaps, the whole sequences are matched.    , global_score, global_align +   -- * Overlapping alignment.+   -- | The suffix of one sequence matches a prefix of another.+   , overlap_score, overlap_align+    -- * Matrix construction    , qualMx++   -- * Interactive testing of alignments+   , test    ) where  import Data.List (maximumBy,partition,unfoldr,zip4,tails) import qualified Data.ByteString.Lazy as B+import qualified Data.ByteString.Lazy.Char8 as BC import Data.Array.Unboxed  import Bio.Sequence.SeqData hiding ((!))@@ -113,13 +125,29 @@ local_score mx g s1 s2 = maximum . map (uncurry max) . concat                          $ columns (score_select 0 mx g) (0,fst g) s1 s2 +-- | Calucalte best overlap score, where gaps at the edges are free+--   The starting point is like for local score (0 cost for initial indels),+--   the result is the maximum anywhere in the last column or bottom row of the matrix.++-- Oh. local score_select will not work, since we cannot replace any matrix entry+-- with zero in order to initiate a new alignment.  So we need 'minf', except in the+-- initial row/column.  Damn.+overlap_score :: QualMx Double -> (Double,Double) -> Sequence -> Sequence -> Double+overlap_score mx g s1 s2 = maximum $ map (uncurry max) $ sel cols+    where cols   = columns overlap_score_select (0,fst g) s1 s2 +          sel cs = map last cs ++ last cs+          -- well - edges have less than three options, so we can set them to zero+          overlap_score_select cds@[_,_,_] = score_select minf mx g cds+          overlap_score_select [_] = (0,minf)+ -- | Generic scoring and selection function for global and local scoring score_select :: Double -> QualMx Double -> (Double,Double) -> QSelector (Double,Double) score_select minf mx (go,ge) cds =     let (reps,ids) = partition (isRepl.snd') cds         s = maximum $ minf:[max sub gap + mx q1 q2 (x,y) | ((sub,gap),Repl x y,q1,q2) <- reps]         g = maximum $ minf:[max (sub+go) (gap+ge) | ((sub,gap),_op,_q1,_q2) <- ids]-    in (s,g)+    in s `seq` g `seq` (s,g) -- (s,g) +             -- seq makes local slower, but overlap faster(!?)  -- ------------------------------------------------------------ -- Alignments (rip from AAlign)@@ -137,19 +165,33 @@ global_align mx g s1 s2 = revsnd . uncurry max' . last . last                $ columns (align_select minf mx g) ((0,[]),(fst g,[])) s1 s2 --- | Calculate local alignmnet (Smith-Waterman)+-- | Calculate local alignment (Smith-Waterman)+--   (can we replace uncurry max' with fst - a local alignment must always end on a subst, no?) local_align :: QualMx Double -> (Double,Double) -> Sequence -> Sequence -> (Double, EditList) local_align mx g s1 s2 = revsnd . maximumBy (compare `on` fst)                          . map (uncurry max') . concat                          $ columns (align_select 0 mx g) ((0,[]),(fst g,[])) s1 s2 -{---- | Calculate the optimal match between a suffix of one sequence and a prefix---   of the other (useful for e.g. sequence assembly)-overlap_align mx g s1 s2 =-     columns (align_select minf mx g) ((0,[]),(fst g,[])) s1 s2--}+-- | Calucalte best overlap score, where gaps at the edges are free+--   The starting point is like for local score (0 cost for initial indels),+--   the result is the maximum anywhere in the last column or bottom row of the matrix.+overlap_align :: QualMx Double -> (Double,Double) -> Sequence -> Sequence -> (Double,EditList)+overlap_align mx g s1 s2 = revsnd . maximumBy (compare `on` fst) . map (uncurry max') $ sel cols+    where cols   = columns overlap_align_select ((0,[]),(minf,[])) s1 s2 +          sel cs = map last cs ++ last cs+          -- again, the edges have less than three options, so we can set them to zero+          overlap_align_select cds@[_,_,_] = align_select minf mx g cds+          overlap_align_select [_] = ((0,[]),(minf,[])) +-- | Variant that retains indels to retain the entire sequence in the result+overlap_align' :: QualMx Double -> (Double,Double) -> Sequence -> Sequence -> (Double,EditList)+overlap_align' mx g s1 s2 = revsnd . maximumBy (compare `on` fst) . map (uncurry max') $ sel cols+    where cols   = columns overlap_align_select ((0,[]),(fst g,[])) s1 s2 +          sel cs = map last cs ++ last cs+          -- the old "fewer choices" trick+          overlap_align_select cds@[_,_,_] = align_select minf mx g cds+          overlap_align_select [(((s1,es1),(s2,es2)),e,_,_)] = ((0,e:es2),(0,e:es2))+ -- (maybe better to reverse the inputs for global?) revsnd (s,a) = (s,reverse a) @@ -164,3 +206,25 @@     in (s,g)  snd' (_,x,_,_) = x++test :: IO ()+test = do+  putStrLn "Enter two strings:"+  s1' <- getLine+  s2' <- getLine+  let s1 = Seq (BC.pack "foo") (BC.pack s1') Nothing+      s2 = Seq (BC.pack "bar") (BC.pack s2') Nothing+      mx = qualMx+      g  = (-5,-2)+  let ga = global_align mx g s1 s2+      la = local_align mx g s1 s2+      oa = overlap_align mx g s1 s2+      or = overlap_align mx g (revcompl s1) (revcompl s2)+  putStrLn ("GLOBAL:   " ++ show (fst ga))+  putStrLn $ showalign $ snd ga+  putStrLn ("OVERLAP:  " ++ show (fst oa))+  putStrLn $ showalign $ snd oa+  putStrLn ("OVERLAP (rc):  " ++ show (fst or))+  putStrLn $ showalign $ snd or+  putStrLn ("LOCAL:    " ++ show (fst la))+  putStrLn $ showalign $ snd la

@@ -10,6 +10,7 @@       Sequence(..), Offset, SeqData, Qual, QualData     -- ** Accessor functions     , seqlength, seqlabel, seqheader, seqdata, seqqual, (!)+    , appendHeader, setHeader      -- ** Converting to and from String.     , fromStr, toStr

@@ -93,9 +93,9 @@   -- warning: assumes correct qual file!   let mkseq s1@(Seq x y _) s2@(Seq _ _ (Just z))           | seqlabel s1 /= seqlabel s2 = error ("mismatching sequence and quality: "-                                                ++show (seqlabel s1,seqlabel s2))+                                                ++toStr (seqlabel s1)++","++toStr (seqlabel s2))           | B.length y /= B.length z   = error ("mismatching sequence and quality lengths:"-                                                ++ show (seqlabel s1,B.length y,B.length z))+                                                ++ toStr (seqlabel s1)++","++show (B.length y,B.length z))           | otherwise   = Seq x y (Just z)       mkseq _ _ = error "readFastaQual: could not combine Fasta and Qual information"   return $ zipWith mkseq ss qs

@@ -1,89 +1,8 @@-{- | -   GOA - parse and index Gene Onthology Annotations-   In particular, the file 'gene_association.goa_uniprot' that contains-   links between GO terms and UniProt accessions.--   (Where to find the hierarchical relationship between GO terms?)-   <http://www.geneontology.org/ontology/gene_ontology.obo> contains isA relationships-   <http://www.geneontology.org/GO.format.obo-1_2.shtml> describes the format--}--module Bio.Sequence.GOA where-import Data.ByteString.Lazy.Char8 (ByteString,pack,unpack,copy)-import qualified Data.ByteString.Lazy.Char8 as B+-- | Moved to GeneOnthology - this is for backwards compatibility. --- | Read the goa_uniprot file (warning: this one is huge!)-readGOA :: FilePath -> IO [Annotation]-readGOA f = B.readFile f >>= -            return . map mkAnn . decomment+module Bio.Sequence.GOA {-# DEPRECATED "Moved to Bio.Sequence.GeneOntology, use that instead." #-}(module GO,readGO) where+import Bio.Sequence.GeneOntology as GO --- | Read GO term definitions readGO :: FilePath -> IO [GoDef]-readGO f = B.readFile f >>= -           return . map mkGoDef . decomment--decomment :: ByteString -> [ByteString]-decomment = filter (\l -> not (B.null l) && B.head l /= '!') . B.lines--newtype GoTerm = GO Int deriving (Eq,Ord)-type UniProtAcc = ByteString-data GoClass = Func | Proc | Comp--instance Read GoTerm where -    readsPrec n ('G':'O':':':xs) = map (\(i,s)-> (GO i,s)) (readsPrec n xs)-    readsPrec n e = error ("couldn't parse GO term: "++show e)-instance Show GoTerm where show (GO x) = "GO:"++show x--instance Read GoClass where-    readsPrec _ ('F':xs) = [(Func,xs)]-    readsPrec _ ('P':xs) = [(Proc,xs)]-    readsPrec _ ('C':xs) = [(Comp,xs)]-    readsPrec _ _ = []-instance Show GoClass where -    show Func = "F"-    show Proc = "P"-    show Comp = "C"---- | GOA Annotation - or multiple annotations?-data Annotation = Ann !UniProtAcc !GoTerm !EvidenceCode deriving (Show)--mkAnn :: ByteString -> Annotation-mkAnn = pick . B.words-    where pick (_db:up:rest) = pick' up $ getGo rest-          pick' up' (go:_:ev:_) = Ann (copy up') (read $ unpack go) (read $ unpack ev)-          getGo = dropWhile (not . B.isPrefixOf (pack "GO:"))---- | GO maps GO terms (GO:xxxx for some number xxxx) to biologically meaningful terms.---   Defined in <http://www.geneontology.org/doc/GO.terms_and_ids>---   The format is  "GO:0000000 [tab] text string  [tab] F|P|C"-data GoDef = GoDef !GoTerm !ByteString !GoClass deriving (Show)--mkGoDef :: ByteString -> GoDef-mkGoDef = pick . B.split '\t'-    where pick [go,desc,cls] = GoDef (read $ unpack go) (copy desc) (read $ unpack cls)-          pick _xs = error ("Couldn't decipher GO definition from: "++show _xs)---- | Evidence codes describe the type of support for an annotation--- <http://www.geneontology.org/GO.evidence.shtml>-data EvidenceCode = IC  -- Inferred by Curator-	          | IDA -- Inferred from Direct Assay-	          | IEA -- Inferred from Electronic Annotation-	          | IEP -- Inferred from Expression Pattern-	          | IGC -- Inferred from Genomic Context-	          | IGI -- Inferred from Genetic Interaction-	          | IMP -- Inferred from Mutant Phenotype-	          | IPI -- Inferred from Physical Interaction-	          | ISS -- Inferred from Sequence or Structural Similarity-	          | NAS -- Non-traceable Author Statement-	          | ND  -- No biological Data available-	          | RCA -- inferred from Reviewed Computational Analysis-	          | TAS -- Traceable Author Statement-	          | NR  -- Not Recorded -     deriving (Read,Show,Eq)---- | The vast majority of GOA data is IEA, while the most reliable information---   is manually curated.  Filtering on this is useful to keep data set sizes---   manageable, too.-isCurated :: EvidenceCode -> Bool-isCurated = not . (`elem` [ND,IEA])+readGO = GO.readTerms 

@@ -0,0 +1,208 @@+{- | +   GeneOntology - parse and index Gene Ontology Annotations+   In particular, the file 'gene_association.goa_uniprot' that contains+   links between GO terms and UniProt accessions.++   * <http://www.geneontology.org/ontology/gene_ontology.obo> +   -- Contains the hierarchy including isA relationships.++   * <http://www.geneontology.org/GO.format.obo-1_2.shtml> +   -- Describes the OBO format.++   * <ftp://ftp.ebi.ac.uk/pub/databases/GO/goa/UNIPROT/> +   -- Contains the GOA-UniProt mapping (and a README file).++   * <http://www.geneontology.org/ontology/GO.defs> +   -- Contains GO definitions (not supported here yet).++   * <http://www.geneontology.org/doc/GO.terms_and_ids> +   -- GO definitions, simpler and more schematically.++-}++module Bio.Sequence.GeneOntology +    (+    -- * Basic data types+    GoTerm(..), GoDef(..)++    -- * Reading the OBO format+    , GoHierarchy, readObo++    -- * Reading 'terms and ids'+    , readTerms++    -- * Reading UniProt associations+    , Annotation(..), UniProtAcc, GoClass(..), EvidenceCode(..), readGOA, isCurated++    -- * Utility stuff+    , decomment++    ) where++import Data.ByteString.Lazy.Char8 (ByteString,pack,unpack,copy)+import qualified Data.ByteString.Lazy.Char8 as B++-- | Read the GO hierarchy from the obo file.  Note that this is not quite a tree structure.+readObo :: FilePath -> IO GoHierarchy+readObo f = B.readFile f >>=+            return . mkGoHier . decomment++-- | Read the goa_uniprot file (warning: this one is huge!)+readGOA :: FilePath -> IO [Annotation]+readGOA f = B.readFile f >>= +            return . map mkAnn . decomment++-- | Read GO term definitions, from the GO.terms_and_ids file+readTerms :: FilePath -> IO [GoDef]+readTerms f = B.readFile f >>= +           return . map mkGoDef . decomment++decomment :: ByteString -> [ByteString]+decomment = filter (\l -> not (B.null l) && B.head l /= '!') . B.lines++-- ----------------------------------------------------------+-- Reading the Obo file containing the ontology definition+-- ----------------------------------------------------------++-- | A list of Go definitions, with pointers to parent nodes.  Read from the .obo file.+--   The user may construct the explicit hierachy by storing these in a Map or similar+type GoHierarchy = [(GoDef,[GoTerm])]++-- Each entry may span multiple lines, thus this function is slightly different from its siblings.+-- Todo: strictness? copy?+mkGoHier :: [ByteString] -> [(GoDef,[GoTerm])]+mkGoHier ls = go $ dropWhile (not . termStart) ls+    where termStart = (== B.pack "[Term]")+          go [] = []+          go (_:zs) = let (this,rest) = span (not . B.isPrefixOf (B.pack "[")) zs in+              if null this then +                  if not (null rest) then error "Parse failure in mkGoHier/go" +                  else []+              else (mk1 $ map ($ this) [getId, getName, getNamespace, getIsA]) : mkGoHier rest++          mk1 xs@[i,n,ns,isa] +              | or (map null [i,n,ns]) = error ("Failed to parse Go Term (missing field in entry):\n"+                                                ++unlines (map unpack $ concat xs))+              | length i /= 1 || length n /= 1 || length ns /= 1 +                         = error ("Failed to parse Go Term (incorrect field multiplicity):\n"+                                  ++unlines (map unpack $ concat xs))+              | otherwise = (GoDef (getGo $ head i) (head n) (readNS $ head ns), map getGo isa)+          mk1 _ = error "This shouldn't happen!"++          getId        = map ((!!1) . B.words) . filter (B.isPrefixOf (pack "id:"))+          getName      = map (B.unwords. tail  . B.words) . filter (B.isPrefixOf (pack "name:"))+          getNamespace = map ((!!1) . B.words) . filter (B.isPrefixOf (pack "namespace:"))+          getIsA       = map ((!!1) . B.words) . filter (B.isPrefixOf (pack "is_a:"))+          readNS xs = case unpack xs of "biological_process" -> Proc+                                        "molecular_function" -> Func+                                        "cellular_component" -> Comp+                                        _ -> error ("Unknown function: "++unpack xs)+++-- ----------------------------------------------------------+-- Reading GoTerms from the GO.terms_and_ids file+-- ----------------------------------------------------------++-- | A GO term is a positive integer+newtype GoTerm = GO Int deriving (Eq,Ord)+data GoClass = Func | Proc | Comp++instance Read GoTerm where +    readsPrec n ('G':'O':':':xs) = map (\(i,s)-> (GO i,s)) (readsPrec n xs)+    readsPrec _ e = error ("couldn't parse GO term: "++show e)+instance Show GoTerm where show (GO x) = "GO:"++show x++getGo :: ByteString -> GoTerm+getGo bs = GO $ fst $ maybe e id (B.readInt $ B.drop 3 bs)+    where e = error ("Unable to parse GO term"++unpack bs)++-- | A GoDef maps a "GoTerm" to a description and a "GoClass".+data GoDef = GoDef !GoTerm !ByteString !GoClass deriving (Show)++--   Defined in <http://www.geneontology.org/doc/GO.terms_and_ids>+--   The format is  "GO:0000000 [tab] text string  [tab] F|P|C"++-- | Parse a "GoDef" from a line in the GO.terms_and_ids file.+mkGoDef :: ByteString -> GoDef+mkGoDef = pick . B.split '\t'+    where pick [go,desc,cls] = GoDef (read $ unpack go) (copy desc) (read $ unpack cls)+          pick _xs = error ("Couldn't decipher GO definition from: "++show _xs)++instance Read GoClass where+    readsPrec _ ('F':xs) = [(Func,xs)]+    readsPrec _ ('P':xs) = [(Proc,xs)]+    readsPrec _ ('C':xs) = [(Comp,xs)]+    readsPrec _ _ = []+instance Show GoClass where +    show Func = "F"+    show Proc = "P"+    show Comp = "C"++-- ----------------------------------------------------------+-- Reading Annotations from the GOA UniProt-GO association file+-- ----------------------------------------------------------++-- | A UniProt identifier (short string of capitals and numbers).+type UniProtAcc = ByteString+-- | A GOA annotation, containing a UniProt identifier, a GoTerm and an evidence code.+data Annotation = Ann !UniProtAcc !GoTerm !EvidenceCode deriving (Show)++-- | Reading an "Annotation" from a line in the association file.+mkAnn :: ByteString -> Annotation+mkAnn = pick . B.words+    where pick (_db:up:rest) = pick' up $ findGo rest+          pick _ = error "Internal error: mkAnn/pick"+          pick' up' (go:_:ev:_) = Ann (copy up') (getGo go) (getEC ev)+          pick' _ _ = error "Internal error: mkAnn/pick'"+          findGo = dropWhile (not . B.isPrefixOf (pack "GO:"))++-- | Evidence codes describe the type of support for an annotation+-- <http://www.geneontology.org/GO.evidence.shtml>+data EvidenceCode = IC  -- ^ Inferred by Curator+	          | IDA -- ^ Inferred from Direct Assay+	          | IEA -- ^ Inferred from Electronic Annotation+	          | IEP -- ^ Inferred from Expression Pattern+	          | IGC -- ^ Inferred from Genomic Context+	          | IGI -- ^ Inferred from Genetic Interaction+	          | IMP -- ^ Inferred from Mutant Phenotype+	          | IPI -- ^ Inferred from Physical Interaction+	          | ISS -- ^ Inferred from Sequence or Structural Similarity+	          | NAS -- ^ Non-traceable Author Statement+	          | ND  -- ^ No biological Data available+	          | RCA -- ^ Inferred from Reviewed Computational Analysis+	          | TAS -- ^ Traceable Author Statement+	          | NR  -- ^ Not Recorded +     deriving (Read,Show,Eq)++-- | Read the evidence code from a ByteString (no error checking!).+getEC :: ByteString -> EvidenceCode +getEC s = case B.uncons s of+            Just ('I',s') -> case B.uncons s' of+                               Just ('C',_) -> IC+                               Just ('D',_) -> IDA+                               Just ('E',s'') -> case B.head s'' of 'A' -> IEA+                                                                    'P' -> IEP+                                                                    _ -> e 1+                               Just ('G',s'') -> case B.head s'' of 'C' -> IGC+                                                                    'I' -> IGI+                                                                    _ -> e 2+                               Just ('M',_) -> IMP+                               Just ('P',_) -> IPI+                               Just ('S',_) -> ISS+                               _ -> e 3+            Just ('N',s') -> case B.head s' of 'A' -> NAS+                                               'D' -> ND+                                               'R' -> NR+                                               _ -> e 4+            Just ('R',_) -> RCA+            Just ('T',_) -> TAS+            _ -> e 5+    where e :: Int -> a+          e n = error ("Illegal GO evidence code ("++show n++"): "++unpack s)++-- | The vast majority of GOA data is IEA, while the most reliable information+--   is manually curated.  Filtering on this is useful to keep data set sizes+--   manageable, too.+isCurated :: EvidenceCode -> Bool+isCurated = not . (`elem` [ND,IEA])+

@@ -0,0 +1,18 @@+{- | +   Functionality for manipulating KEGG annotations.++   KEGG is a bit hard find, but there exist species-specific tables+   Available organisms are listed in the table at++      <ftp://ftp.genome.jp/pub/kegg/genes/etc/all_species.tab>++   Data for each organism is stored its own subdirectory under++      <ftp://ftp.genome.jp/pub/kegg/genes/organisms/>++   Containing tables linking everything -- including external resources like+   UniProt, PDB, or NCBI -- together.+-}++module Bio.Sequence.KEGG where+

@@ -19,6 +19,9 @@     , (!), seqlength,seqlabel,seqheader,seqdata     , (?), hasqual, seqqual +      -- * Adding information to header+    , appendHeader, setHeader+       -- * Converting to and from [Char]     , fromStr, toStr @@ -75,6 +78,8 @@ {-# INLINE (?) #-} (?) :: Sequence -> Offset -> Qual (?) (Seq _ _ (Just qs)) = BB.index qs+(?) (Seq _ _ Nothing) = +    error "Attempting to use 'seqqual' on sequence without associated quality data"  -- | Return sequence length. seqlength :: Sequence -> Offset@@ -95,10 +100,18 @@ -- | Return the quality data, or error if none exist.  Use hasqual if in doubt. seqqual :: Sequence -> QualData seqqual (Seq _ _ (Just q)) = q+seqqual (Seq _ _ Nothing) = +    error "Attempting to use 'seqqual' on sequence without associated quality data"  -- | Check whether the sequence has associated quality data. hasqual :: Sequence -> Bool hasqual (Seq _ _ q) = isJust q++-- | Modify the header by appending text, or by replacing+--   all but the sequence label (i.e. first word).+appendHeader, setHeader :: Sequence -> String -> Sequence+appendHeader (Seq h d q) t = (Seq (B.append h (B.pack (" "++t))) d q)+setHeader (Seq h d q) t = (Seq (B.append (head $ B.words h) (B.pack (" "++t))) d q)  ------------------------------------------------------------ -- Nucleotide (DNA, RNA) specific stuff

@@ -2,27 +2,28 @@    Utility module, with various useful stuff. -} -module Bio.Util where+module Bio.Util (lines, splitWhen, countIO, sequence', mylines) where +import Prelude hiding (lines) import System.IO        (stderr, hPutStr, hFlush) import System.IO.Unsafe (unsafeInterleaveIO) import Data.List        (groupBy) import qualified Data.ByteString.Lazy.Char8 as B --- workaround, current ByteString.Lazy.Char8 contains a bug in 'lines'+-- | Workaround, the current "Data.ByteString.Lazy.Char8" contains a bug in 'Data.ByteString.Lazy.Char8.lines'.+lines, mylines :: B.ByteString -> [B.ByteString] lines = case length (B.lines $ B.pack "\n") of 1 -> B.lines                                                0 -> mylines mylines s = case B.elemIndex '\n' s of               Nothing -> if B.null s then [] else [s]               Just i  -> B.take i s : mylines (B.drop (i+1) s) --- | Break a list of bytestrings on a predicate+-- | Break a list of bytestrings on a predicate. splitWhen :: (B.ByteString -> Bool) -> [B.ByteString] -> [[B.ByteString]] splitWhen p = groupBy (\_ y -> B.null y || not (p y)) --- | Output (to stderr) progress while evaluating a lazy list+-- | Output (to stderr) progress while evaluating a lazy list. --   Useful for generating output while (conceptually, at least) in pure code---   Strictness warning!! This doesn't *quite* work in all cases.  Why? countIO :: String -> String -> Int -> [a] -> IO [a] countIO msg post step xs = sequence' $ map unsafeInterleaveIO ((blank >> outmsg (0::Int) >> c):cs)    where (c:cs) = ct 0 xs@@ -35,7 +36,7 @@                                 []   -> map return (init a) ++ [outmsg (s+length a) >> hPutStr stderr post >> return (last a)]                                 _ -> map return a ++ [outmsg s >> return (head b)] ++ ct next (tail b) --- lazy sequence+-- | A lazier version of 'Control.Monad.sequence' in "Control.Monad", needed by 'countIO' above. sequence' :: [IO a] -> IO [a] sequence' ms = foldr k (return []) ms     where k m m' = do { x <- m; xs <- unsafeInterleaveIO m'; return (x:xs) }

@@ -1,4 +1,4 @@--- | Lazy "many" combinator for Parsec.+-- | Lazy \"many\" combinator for Parsec. --   Courtesy of Tomasz Zielonka.  module Bio.Util.Parsex where

@@ -1,5 +1,5 @@ Name:                bio-Version:             0.3.3.2+Version:             0.3.3.4 License:             LGPL License-file:        LICENSE Author:              Ketil Malde@@ -36,6 +36,8 @@                      Bio.Sequence.Fasta, Bio.Sequence.TwoBit, Bio.Sequence.Phd,                      Bio.Sequence.Entropy, Bio.Sequence.HashWord,                      Bio.Sequence.GOA,+                     Bio.Sequence.GeneOntology,+		     Bio.Sequence.KEGG,                      Bio.Alignment.BlastData, Bio.Alignment.BlastFlat,                      Bio.Alignment.Blast, Bio.Alignment.BlastXML,                      Bio.Alignment.AlignData, Bio.Alignment.Matrices,