diff --git a/Bio/Alignment/ACE.hs b/Bio/Alignment/ACE.hs --- a/Bio/Alignment/ACE.hs +++ b/Bio/Alignment/ACE.hs @@ -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 diff --git a/Bio/Alignment/Blast.hs b/Bio/Alignment/Blast.hs --- a/Bio/Alignment/Blast.hs +++ b/Bio/Alignment/Blast.hs @@ -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 "" -- 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 diff --git a/Bio/Alignment/BlastFlat.hs b/Bio/Alignment/BlastFlat.hs --- a/Bio/Alignment/BlastFlat.hs +++ b/Bio/Alignment/BlastFlat.hs @@ -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.) diff --git a/Bio/Alignment/BlastXML.hs b/Bio/Alignment/BlastXML.hs --- a/Bio/Alignment/BlastXML.hs +++ b/Bio/Alignment/BlastXML.hs @@ -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 diff --git a/Bio/Alignment/Matrices.hs b/Bio/Alignment/Matrices.hs --- a/Bio/Alignment/Matrices.hs +++ b/Bio/Alignment/Matrices.hs @@ -4,28 +4,39 @@ When in doubt, use BLOSUM62. Consult for some hints on good parameters for nucleotide alignments. + + See also 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), diff --git a/Bio/Alignment/QAlign.hs b/Bio/Alignment/QAlign.hs --- a/Bio/Alignment/QAlign.hs +++ b/Bio/Alignment/QAlign.hs @@ -8,6 +8,8 @@ This module performs sequences alignments, takes quality values into account. + + See also . -} {-# 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 diff --git a/Bio/Sequence.hs b/Bio/Sequence.hs --- a/Bio/Sequence.hs +++ b/Bio/Sequence.hs @@ -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 diff --git a/Bio/Sequence/Fasta.hs b/Bio/Sequence/Fasta.hs --- a/Bio/Sequence/Fasta.hs +++ b/Bio/Sequence/Fasta.hs @@ -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 diff --git a/Bio/Sequence/GOA.hs b/Bio/Sequence/GOA.hs --- a/Bio/Sequence/GOA.hs +++ b/Bio/Sequence/GOA.hs @@ -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?) - contains isA relationships - 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 --- 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 --- -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 diff --git a/Bio/Sequence/GeneOntology.hs b/Bio/Sequence/GeneOntology.hs new file mode 100644 --- /dev/null +++ b/Bio/Sequence/GeneOntology.hs @@ -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. + + * + -- Contains the hierarchy including isA relationships. + + * + -- Describes the OBO format. + + * + -- Contains the GOA-UniProt mapping (and a README file). + + * + -- Contains GO definitions (not supported here yet). + + * + -- 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 +-- 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 +-- +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]) + diff --git a/Bio/Sequence/KEGG.hs b/Bio/Sequence/KEGG.hs new file mode 100644 --- /dev/null +++ b/Bio/Sequence/KEGG.hs @@ -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 + + + + Data for each organism is stored its own subdirectory under + + + + Containing tables linking everything -- including external resources like + UniProt, PDB, or NCBI -- together. +-} + +module Bio.Sequence.KEGG where + diff --git a/Bio/Sequence/SeqData.hs b/Bio/Sequence/SeqData.hs --- a/Bio/Sequence/SeqData.hs +++ b/Bio/Sequence/SeqData.hs @@ -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 diff --git a/Bio/Util.hs b/Bio/Util.hs --- a/Bio/Util.hs +++ b/Bio/Util.hs @@ -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) } diff --git a/Bio/Util/Parsex.hs b/Bio/Util/Parsex.hs --- a/Bio/Util/Parsex.hs +++ b/Bio/Util/Parsex.hs @@ -1,4 +1,4 @@ --- | Lazy "many" combinator for Parsec. +-- | Lazy \"many\" combinator for Parsec. -- Courtesy of Tomasz Zielonka. module Bio.Util.Parsex where diff --git a/bio.cabal b/bio.cabal --- a/bio.cabal +++ b/bio.cabal @@ -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,