bio 0.3.3.4 → 0.3.5

raw patch · 26 files changed

+1731/−61 lines, 26 filesdep +mtldep +old-timedep +randomdep ~bytestringPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: mtl, old-time, random

Dependency ranges changed: bytestring

API changes (from Hackage documentation)

Files

• Bio/Alignment/Blast.hs +2/−0
• Bio/Alignment/BlastXML.hs +8/−2
• Bio/Alignment/Multiple.hs +8/−8
• Bio/Alignment/Soap.hs +134/−0
• Bio/GFF3/Escape.hs +55/−0
• Bio/GFF3/Feature.hs +139/−0
• Bio/GFF3/FeatureHier.hs +115/−0
• Bio/GFF3/FeatureHierSequences.hs +87/−0
• Bio/GFF3/SGD.hs +94/−0
• Bio/Location/ContigLocation.hs +135/−0
• Bio/Location/LocMap.hs +115/−0
• Bio/Location/Location.hs +119/−0
• Bio/Location/OnSeq.hs +48/−0
• Bio/Location/Position.hs +41/−0
• Bio/Location/SeqLocMap.hs +33/−0
• Bio/Location/SeqLocation.hs +46/−0
• Bio/Location/Strand.hs +33/−0
• Bio/Sequence.hs +6/−0
• Bio/Sequence/FastQ.hs +70/−0
• Bio/Sequence/Fasta.hs +9/−38
• Bio/Sequence/KEGG.hs +30/−0
• Bio/Sequence/SFF.hs +247/−0
• Bio/Sequence/TwoBit.hs +5/−4
• Bio/Util.hs +2/−0
• Bio/Util/TestBase.hs +132/−0
• bio.cabal +18/−9

@@ -73,6 +73,7 @@  -- parse metadata from the preamble parse_preamble :: [ByteString] -> BlastResult+parse_preamble [] = error "Blast: can't parse preamble - empty input" 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@@ -108,3 +109,4 @@                  q_from = ee, q_to = ee, h_from = ee, h_to = ee                }      where ee = error "hit coordinates is not supported - use XML format"+parse_match _ = error "Blast: parse_match: pattern match failed"

@@ -18,12 +18,17 @@ import Text.HTML.TagSoup import Control.Monad import Control.Parallel+import Data.List (isPrefixOf)  -- | Parse BLAST results in XML format readXML :: FilePath -> IO [BlastResult] readXML fp = do -    ts <- return . parseTags =<< readFile fp-    let (h:iters) = breaks (\t -> isTagOpenName "Iteration" t || isTagOpenName "Hit" t) ts+    fc <- readFile fp+    when (not $ isPrefixOf "<?xml" fc) +             $ error ("Bio.Sequence.BlastXML.readXML:\n   The file '"+                      ++fp++"' does not look like an XML file - aborting!")+    let ts = parseTags fc+        (h:iters) = breaks (\t -> isTagOpenName "Iteration" t || isTagOpenName "Hit" t) ts     return [xml2br h iters]  -- | breaks p = groupBy (const (not.p))@@ -98,6 +103,7 @@     where get = getFrom ms           mkFrame f = Frame (strand' $ signum $ read f) (abs $ read f)           mkStrands h q = Strands (strand' $ read h) (strand' $ read q)+          strand' :: Int -> Strand           strand' s = case s of 1 -> Plus; -1 -> Minus                                 _ -> error ("Strand must be +1 or -1, but was: "++show s) 

@@ -9,7 +9,7 @@  import Bio.Alignment.AlignData import Bio.Sequence-import Bio.Clustering+-- import Bio.Clustering  -- | Progressive multiple alignment. --   Calculate a tree from agglomerative clustering, then align@@ -21,15 +21,15 @@ --    This is central for 'Coffee' evaluation of alignments, and T-Coffee construction --    of alignments. indirect :: EditList -> EditList -> EditList-indirect (Repl x1 x2:xs) (Repl y1 y2:ys) = Repl x1 y2 : indirect xs ys -- assert x2==y1-indirect xs@(Repl _ _:_) (Ins y1:ys)     = Ins y1     : indirect xs ys-indirect (Repl x1 _:xs) (Del y1:ys)      = Del x1     : indirect xs ys+indirect (Repl x1 _x2:xs) (Repl _y1 y2:ys) = Repl x1 y2 : indirect xs ys -- assert x2==y1+indirect xs@(Repl _ _:_) (Ins y1:ys)       = Ins y1     : indirect xs ys+indirect (Repl x1 _:xs) (Del _y1:ys)       = Del x1     : indirect xs ys -indirect (Del x1:xs) ys                  = Del x1     : indirect xs ys -- imply del+ins/=repl+indirect (Del x1:xs) ys                    = Del x1     : indirect xs ys -- imply del+ins/=repl -indirect (Ins x1:xs) (Repl _ y2:ys)      = Ins y2     : indirect xs ys -- assert x1==y1-indirect (Ins x1:xs) (Del y1:ys)         = indirect xs ys -- assert x1 == y1-indirect xs@(Ins _:_) (Ins y1:ys)        = Ins y1 : indirect xs ys+indirect (Ins _x1:xs) (Repl _ y2:ys)       = Ins y2     : indirect xs ys -- assert x1==y1+indirect (Ins _x1:xs) (Del _y1:ys)         = indirect xs ys -- assert x1 == y1+indirect xs@(Ins _:_) (Ins y1:ys)          = Ins y1 : indirect xs ys  indirect [] ys                           = ys -- assert: all Ins indirect xs []                           = xs -- assert: all Del

@@ -0,0 +1,134 @@+module Bio.Alignment.Soap ( SoapAlign(..), SoapAlignMismatch(..)+                          , refSeqPos, refCSeqLoc, refSeqLoc, mismatchSeqPos+                          , parse, unparse, parseMismatch, unparseMismatch+                          , group+                          )+    where ++import Prelude hiding (length)+import Control.Monad.Error+import qualified Data.ByteString.Lazy as LBSW+import qualified Data.ByteString.Lazy.Char8 as LBS+import Data.Char+import qualified Data.List as List (length)+import Data.List (groupBy)++import qualified Bio.Location.ContigLocation as CLoc+import qualified Bio.Location.Location as Loc+import qualified Bio.Location.Position as Pos+import Bio.Location.OnSeq+import qualified Bio.Location.SeqLocation as SeqLoc+import Bio.Location.Strand+import Bio.Sequence.SeqData+++-- | Alignment output from SOAP+data SoapAlign = SA { name :: !SeqName+                    , sequ :: !SeqData     -- ^ Reference strand orientation sequence+                    , qual :: !QualData    -- ^ Reference strand orientation quality data+                    , nhit :: !Int+                    , pairend :: !Char+                    , length :: !Offset+                    , strand :: !Strand+                    , refname :: !SeqName+                    , refstart :: !Offset  -- ^ 1-based index, as output by SOAP, of reference strand 5' end+                    , nmismatch :: !Int+                    , mismatches :: ![SoapAlignMismatch]+                    } deriving (Read, Show, Eq, Ord)++data SoapAlignMismatch = SAM { readnt :: !Char   -- ^ Read nt in reference strand orientation+                             , refnt :: !Char    -- ^ Reference nt in reference strand orientation+                             , offset :: !Offset -- ^ Offset from reference strand 5' end in reference strand orientation+                             , qualnt :: !Qual   -- ^ Quality score of read nt+                             } deriving (Read, Show, Eq, Ord)++mismatchSeqPos :: SoapAlign -> SoapAlignMismatch -> SeqLoc.SeqPos+mismatchSeqPos sa sam = OnSeq (refname sa) $ Pos.Pos (refstart sa + offset sam - 1) (strand sa)++refSeqPos :: SoapAlign -> SeqLoc.SeqPos+refSeqPos sa = OnSeq (refname sa) $ CLoc.startPos $ refCLoc sa++refCSeqLoc :: SoapAlign -> SeqLoc.ContigSeqLoc+refCSeqLoc sa = OnSeq (refname sa) (refCLoc sa)++refCLoc :: SoapAlign -> CLoc.ContigLoc+refCLoc sa = CLoc.ContigLoc (refstart sa - 1) (length sa) (strand sa)++refSeqLoc :: SoapAlign -> SeqLoc.SeqLoc+refSeqLoc sa = OnSeq (refname sa) (Loc.Loc [ refCLoc sa ])++qualScale :: Qual+qualScale = 64++parse :: (Error e, MonadError e m) => LBS.ByteString -> m SoapAlign+parse bstr+    = case LBS.split '\t' bstr of+        (nameStr:sequStr:qualStr:nhitStr:pairendStr:lengthStr:strandStr:refnameStr:refstartStr:nmismatchStr:mismatchesStrs)+            -> do let scQualStr = LBSW.map (subtract qualScale) qualStr+                  nhitInt <- parseIntBStr nhitStr+                  pairendCh <- parseCharBStr pairendStr+                  lengthOff <- parseOffsetBStr lengthStr+                  strandVal <- parseCharBStr strandStr >>= parseStrandChar+                  refstartOff <- parseOffsetBStr refstartStr+                  nmismatchInt <- parseIntBStr nmismatchStr+                  mismatchesList <- mapM parseMismatch mismatchesStrs+                  verifyNMismatches nmismatchInt mismatchesList+                  return $ SA nameStr sequStr scQualStr nhitInt pairendCh lengthOff strandVal refnameStr refstartOff nmismatchInt mismatchesList+        fs -> throwError $ strMsg $ "Bio.Alignment.Soap.parse: too few fields in soapAlign: " ++ show (List.length fs)+    where verifyNMismatches n l = unless (List.length l == n) $ throwError $ strMsg $+                                  "Bio.Alignment.Soap.parse: wrong number of errors: " ++ show (n, List.length l)+          parseStrandChar '+' = return Fwd+          parseStrandChar '-' = return RevCompl+          parseStrandChar ch  = throwError $ strMsg $ "Unknown strand " ++ show ch++parseMismatch :: (Error e, MonadError e m) => LBS.ByteString -> m SoapAlignMismatch+parseMismatch mstr = do (ref, rest1) <- maybe malformed return $ LBS.uncons mstr+                        unless (LBS.isPrefixOf (LBS.pack "->") rest1) malformed+                        let (offstr, rest2) = LBS.span isDigit $ LBS.drop 2 rest1+                        (rd, qualstr) <- maybe malformed return $ LBS.uncons rest2+                        o <- parseOffsetBStr offstr+                        q <- liftM fromIntegral $ parseIntBStr qualstr+                        return $ SAM rd ref o q+    where malformed = throwError $ strMsg $ "Bio.Alignment.Soap.parseMismatch: malformed mismatch field " +                      ++ (show $ LBS.unpack mstr)++parseIntBStr :: (Error e, MonadError e m) => LBS.ByteString -> m Int+parseIntBStr zstr = case LBS.readInt zstr of+                      Just (z, rest) | LBS.null rest -> return z+                      _ -> throwError $ strMsg $ "parseIntBStr: Malformed int field " ++ show zstr++parseOffsetBStr :: (Error e, MonadError e m) => LBS.ByteString -> m Offset+parseOffsetBStr zstr = case LBS.readInteger zstr of+                         Just (z, rest) | LBS.null rest -> return $ fromIntegral z+                         _ -> throwError $ strMsg $ "parseOffsetBStr: Malformed integer field " ++ show zstr++parseCharBStr :: (Error e, MonadError e m) => LBS.ByteString -> m Char+parseCharBStr chstr | LBS.length chstr == 1 = return $ LBS.head chstr+                    | otherwise = throwError $ strMsg $ "parseCharBStr: Malformed char field " ++ show chstr++unparse :: SoapAlign -> LBS.ByteString+unparse sa = LBS.intercalate (LBS.singleton '\t') $ [ name sa+                                                    , sequ sa+                                                    , LBSW.map (+ qualScale) $ qual sa+                                                    , LBS.pack $ show $ nhit sa+                                                    , LBS.singleton $ pairend sa+                                                    , LBS.pack $ show $ length sa+                                                    , unparseStrand $ strand sa+                                                    , refname sa+                                                    , LBS.pack $ show $ refstart sa+                                                    , LBS.pack $ show $ nmismatch sa+                                                    ] ++ map unparseMismatch (mismatches sa)+    where unparseStrand Fwd      = LBS.singleton '+'+          unparseStrand RevCompl = LBS.singleton '-'++unparseMismatch :: SoapAlignMismatch -> LBS.ByteString+unparseMismatch (SAM rd ref off q) = LBS.concat [ LBS.singleton ref+                                                , LBS.pack "->"+                                                , LBS.pack $ show off+                                                , LBS.singleton rd+                                                , LBS.pack $ show q+                                                ]++group :: [SoapAlign] -> [[SoapAlign]]+group = groupBy (equating name)+    where equating f x y = (f x) == (f y)

@@ -0,0 +1,55 @@+module Bio.GFF3.Escape ( unEscapeByteString+                       , escapeByteString, escapeAllBut, escapeAllOf+                       )+    where++import Control.Monad.Error+import Data.Array.ST+import Data.Array.Unboxed+import qualified Data.ByteString.Lazy as LBSW+import qualified Data.ByteString.Lazy.Char8 as LBS+import Data.Char+import Data.Word++-- (Un)Escaping++unEscapeByteString :: (Error e, MonadError e m) => LBS.ByteString -> m LBS.ByteString+unEscapeByteString bstr = case LBS.split '%' bstr of+                            [] -> return LBS.empty+                            [_] -> return bstr+                            (bstr0:brest) -> foldM unEscapeOneSplit bstr0 brest+    where unEscapeOneSplit leftBStr nextBStr+              = do (x1, next1) <- maybe shortEscape return $ LBS.uncons nextBStr+                   (x2, rest) <- maybe shortEscape return $ LBS.uncons next1+                   xch <- unescapeChar x1 x2+                   return $ leftBStr `LBS.append` LBS.singleton xch `LBS.append` rest+          unescapeChar x1 x2+              | isHexDigit x1 && isHexDigit x2 = return $ chr (digitToInt x1 * 16 + digitToInt x2)+              | otherwise = throwError $ strMsg $ "Bad %-escape, " ++ show ['%',x1,x2]+          shortEscape = throwError $ strMsg $ "Bad %-escape, too short"++escapeByteString :: (Char -> Bool) -> LBS.ByteString -> LBS.ByteString+escapeByteString isNotEscaped bstr +    = case LBS.span isNotEscaped bstr of+        (leftStr, rightStr) -> case LBSW.uncons rightStr of+                                 Nothing -> leftStr+                                 Just (escch, rest) -> leftStr `LBS.append` escapeWord8 escch `LBS.append` escapeByteString isNotEscaped rest++escapeWord8 :: Word8 -> LBS.ByteString+escapeWord8 wch = case (fromIntegral wch) `quotRem` 16 of+                    (lch, rch) -> LBS.pack ['%',intToDigit lch, intToDigit rch]+                                 +escapeTable :: Bool -> [Char] -> UArray Char Bool+escapeTable def exc = runSTUArray $ do starr <- newArray (minBound, maxBound) def+                                       mapM_ (\ch -> writeArray starr ch $ not def) exc+                                       return starr++escapeAllBut :: String -> LBS.ByteString -> LBS.ByteString+escapeAllBut notEscaped = escapeByteString isNotEscaped+    where isNotEscaped ch = table ! ch+          table = escapeTable False notEscaped++escapeAllOf :: String -> LBS.ByteString -> LBS.ByteString+escapeAllOf allEscaped = escapeByteString isNotEscaped+    where isNotEscaped ch = table ! ch+          table = escapeTable True allEscaped

@@ -0,0 +1,139 @@+module Bio.GFF3.Feature ( GFFAttr(..)+                        , Feature(..), length+                        , parse, unparse+                        , parseWithFasta+                        , attrByTag+                        , ids, parentIds+                        , contigLoc, loc, seqLoc+                        , name+                        )+    where++import Prelude hiding (length)+import Control.Monad.Error+import qualified Data.ByteString.Lazy.Char8 as LBS+import Data.List hiding (length)+import qualified Data.List as List (length)+import Data.Maybe++import Bio.GFF3.Escape++import qualified Bio.Location.ContigLocation as CLoc+import qualified Bio.Location.Location as Loc+import Bio.Location.OnSeq+import qualified Bio.Location.SeqLocation as SeqLoc+import Bio.Location.Strand+import Bio.Sequence.SeqData++data GFFAttr = GFFAttr { attrTag :: !LBS.ByteString, attrValues :: ![LBS.ByteString] } deriving (Eq, Ord, Show)++data Feature = Feature { seqid :: !LBS.ByteString+                       , source :: !LBS.ByteString+                       , ftype :: !LBS.ByteString+                       , start, end :: !Offset+                       , score :: !(Maybe Double)+                       , strand :: !(Maybe Strand)+                       , phase :: !(Maybe Offset)+                       , attributes :: ![GFFAttr]+                       } deriving (Eq, Ord, Show)++length :: Feature -> Offset+length f = 1 + end f - start f++dot :: LBS.ByteString+dot = LBS.singleton '.' -- for Nothing in optional fields++parse :: (Error e, MonadError e m) => LBS.ByteString -> m Feature+parse line = case LBS.split '\t' line of+               [seqidStr,sourceStr,typeStr,startStr,endStr,scoreStr,strandStr,phaseStr,attrStr]+                   -> do fSeqid <- unEscapeByteString seqidStr+                         fSource <- unEscapeByteString sourceStr+                         fFtype <- unEscapeByteString typeStr+                         fStart <- unEscapeByteString startStr >>= parseInteger+                         fEnd <- unEscapeByteString endStr >>= parseInteger+                         fScore <- unEscapeByteString scoreStr >>= parseMaybe parseDouble+                         fStrand <- unEscapeByteString strandStr >>= parseMaybe parseStrand+                         fPhase <- unEscapeByteString phaseStr >>= parseMaybe parseInteger+                         fAttrs <- parseAttrs attrStr+                         return $ Feature fSeqid fSource fFtype fStart fEnd fScore fStrand fPhase fAttrs+               fs -> throwError $ strMsg $ "Malformed GFF line with " ++ show (List.length fs) ++ " tab-delimited fields"+    where parseInteger zstr = case LBS.readInteger zstr of+                                Just (z, rest) | LBS.null rest -> return $ fromIntegral z+                                _ ->  throwError $ strMsg $ "Malformed integer " ++ show (LBS.unpack zstr)+          parseMaybe subparse str = if str == dot then return Nothing else liftM Just $ subparse str+          parseDouble xstr = case reads $ LBS.unpack xstr of+                               [(x, "")] -> return x+                               _ -> throwError $ strMsg $ "Malformed double " ++ show (LBS.unpack xstr)+          parseStrand sstr | sstr == (LBS.singleton '+') = return Fwd+                           | sstr == (LBS.singleton '-') = return RevCompl+                           | otherwise = throwError $ strMsg $ "Malformed strand " ++ show (LBS.unpack sstr)+          parseAttrs = mapM parseAttr . LBS.split ';'+          parseAttr attrStr = case LBS.break (== '=') attrStr of+                                (tagStr,equalValStr) +                                    -> case LBS.uncons equalValStr of+                                         Just ('=', valStr) -> do tag <- unEscapeByteString tagStr+                                                                  vals <- parseAttrVals valStr+                                                                  return $ GFFAttr tag vals+                                         _ -> throwError $ strMsg $ "Malformed attribute " ++ show (LBS.unpack attrStr)+          parseAttrVals valStr | LBS.null valStr = return [LBS.empty] -- Special case for single null value+                               | otherwise = mapM unEscapeByteString $ LBS.split ',' valStr++escapeField :: LBS.ByteString -> LBS.ByteString+escapeField = escapeAllOf "\t\r\n;=%&,"++escapeSeqid :: LBS.ByteString -> LBS.ByteString+escapeSeqid = escapeAllBut $ ".:^*$@!+_?-|" ++ ['a'..'z'] ++ ['A'..'Z'] ++ ['0'..'9']++unparse :: Feature -> LBS.ByteString+unparse feat = LBS.intercalate (LBS.singleton '\t') [ escapeSeqid $ seqid feat+                                                    , escapeField $ source feat+                                                    , escapeField $ ftype feat+                                                    , escapeField $ unparseInteger $ start feat+                                                    , escapeField $ unparseInteger $ end feat+                                                    , escapeField $ unparseMaybe (LBS.pack . show) $ score feat+                                                    , escapeField $ unparseMaybe unparseStrand $ strand feat+                                                    , escapeField $ unparseMaybe (LBS.pack . show) $ phase feat+                                                    , unparseAttributes $ attributes feat+                                                    ]+    where unparseMaybe unparseJust = maybe dot unparseJust+          unparseStrand Fwd = LBS.singleton '+'+          unparseStrand RevCompl = LBS.singleton '-'+          unparseInteger = LBS.pack . show+          unparseAttributes = LBS.intercalate (LBS.singleton ';') . map unparseAttribute+          unparseAttribute (GFFAttr tag values) = LBS.concat [ escapeField tag, LBS.singleton '=', +                                                               LBS.intercalate (LBS.singleton ',') $ map escapeField values ]++gffFastaDirective :: LBS.ByteString+gffFastaDirective = LBS.pack "##FASTA"++parseWithFasta :: (Error e, MonadError e m) => LBS.ByteString -> m ([Feature], [LBS.ByteString])+parseWithFasta str = case break (== gffFastaDirective) $ LBS.lines str of+                       (featStrs, rest) -> do gffLines <- mapM parse $ filter notComment featStrs+                                              return (gffLines, rest)+    where notComment = maybe False ((/= '#') . fst) . LBS.uncons ++attrByTag :: LBS.ByteString -> Feature -> [LBS.ByteString]+attrByTag tag = findByTag . attributes+    where findByTag = maybe [] attrValues . find ((== tag) . attrTag)++idTag, parentTag :: LBS.ByteString+idTag = LBS.pack "ID"+parentTag = LBS.pack "Parent"++ids :: Feature -> [LBS.ByteString]+ids = attrByTag idTag++parentIds :: Feature -> [LBS.ByteString]+parentIds = attrByTag parentTag++contigLoc :: Feature -> CLoc.ContigLoc+contigLoc f = CLoc.ContigLoc (start f - 1) (1 + end f - start f) (fromMaybe Fwd $ strand f)++loc :: Feature -> Loc.Loc+loc f = Loc.Loc [ contigLoc f ]++seqLoc :: Feature -> SeqLoc.SeqLoc+seqLoc f = OnSeq (seqid f) (loc f)++name :: (Error e, MonadError e m) => Feature -> m SeqName+name f = maybe (throwError $ strMsg $ "No ID for feature " ++ show f) return $ listToMaybe $ ids f

@@ -0,0 +1,115 @@+{-# OPTIONS_GHC -XFlexibleContexts #-}++module Bio.GFF3.FeatureHier ( FeatureHier, features, lookupId, lookupIdChildren+                            , fromList, insert, delete+                            , parents, children, parentsM, childrenM+                            , checkInvariants+                            )+    where ++import Control.Monad.Error+import Control.Monad.Reader+import qualified Data.ByteString.Lazy.Char8 as LBS+import Data.Graph+import Data.List hiding (insert, delete)+import qualified Data.List as List (delete)+import qualified Data.Map as M+import Data.Maybe+import qualified Data.Set as S++import Bio.GFF3.Feature++data FeatureHier = FeatureHier { features :: !(S.Set Feature)+                               , idToFeature :: !(M.Map LBS.ByteString Feature)+                               , idToChildren :: !(M.Map LBS.ByteString [Feature])+                               } deriving (Show)++empty :: FeatureHier+empty = FeatureHier { features = S.empty, idToFeature = M.empty, idToChildren = M.empty }++lookupId :: (Error e, MonadError e m) => FeatureHier -> LBS.ByteString -> m Feature+lookupId fh idstr = maybe idNotFound return $ M.lookup idstr $ idToFeature fh+    where idNotFound = throwError $ strMsg $ "lookupId: feature ID " ++ show (LBS.unpack idstr) ++ " not found"++lookupIdChildren :: (Error e, MonadError e m) => FeatureHier -> LBS.ByteString -> m [Feature]+lookupIdChildren fh idstr = maybe idNotFound return $ M.lookup idstr $ idToChildren fh+    where idNotFound = throwError $ strMsg $ "lookupIdChildren: feature ID " ++ show (LBS.unpack idstr) ++ " not found"++fromList :: (Error e, MonadError e m) => [Feature] -> m FeatureHier+fromList feats = foldM (flip insert) empty $ parentsFirst feats++-- topSort, i < j when i --> j and j /-> i+-- parents first mean edges go parent --> child+parentsFirst :: [Feature] -> [Feature]+parentsFirst feats = let !(graph, fromVertex, _) = graphFromEdges $ featureEdges feats+                     in map (toFeature fromVertex) $ topSort graph+    where toFeature fromVertex v = case fromVertex v of (feat, _, _) -> feat++-- Feature -> (Feature, its key, keys of its out-edges)+featureEdges :: [Feature] -> [(Feature, Int, [Int])]+featureEdges feats = map featureEdge $ keyedFeats+    where keyedFeats = zip feats [1..]+          idToChildKeys = foldl' insertParentsToKey M.empty keyedFeats+          insertParentsToKey m0 (f, k) = foldl' (insertParentToKey k) m0 $ parentIds f+          insertParentToKey k m0 pid = M.insertWith' (++) pid [k] m0+          featureEdge (feat, key) = (feat, key, concatMap childKeys $ ids feat)+          childKeys pid = M.findWithDefault [] pid idToChildKeys+++insert :: (Error e, MonadError e m) => Feature -> FeatureHier -> m FeatureHier+insert f hier0+    = liftM3 FeatureHier doNewAllFeatures doNewIdToFeature doNewIdToChildren+    where doNewAllFeatures = return $ S.insert f $ features hier0+          doNewIdToFeature = foldM insertIdToFeature (idToFeature hier0) $ ids f+          insertIdToFeature m0 fid = if M.member fid m0+                                     then throwError $ strMsg $ "insertFeature: Duplicate ID " ++ show fid+                                     else return $ M.insert fid f m0+          doNewIdToChildren = foldM insertIdToChildren (idToChildren hier0) $ parentIds f+          insertIdToChildren m0 pid = if M.member pid $ idToFeature hier0+                                      then return $ M.insertWith' (++) pid [f] m0+                                      else throwError $ strMsg $ "insertFeature: Parent ID " ++ show pid ++ " not present"++delete :: (Error e, MonadError e m) => Feature -> FeatureHier -> m FeatureHier+delete f hier0+    = liftM3 FeatureHier doNewAllFeatures doNewIdToFeature doNewIdToChildren+    where doNewAllFeatures = if S.member f $ features hier0+                             then return $ S.delete f $ features hier0+                             else throwError $ strMsg $ "deleteFeature: Feature not present " ++ show f+          doNewIdToFeature = foldM deleteIdToFeature (idToFeature hier0) $ ids f+          deleteIdToFeature m0 fid = if M.member fid m0+                                     then return $ M.delete fid m0+                                     else throwError $ strMsg $ "deleteFeature: ID not present for feature " ++ show (fid, f)+          doNewIdToChildren = if any (flip M.member $ idToChildren hier0) $ ids f+                              then throwError $ strMsg $ "deleteFeature: Feature has children: " ++ show f+                              else foldM deleteIdToChildren (idToChildren hier0) $ parentIds f+          deleteIdToChildren m0 pid = if maybe False (elem f) $ M.lookup pid m0+                                      then throwError $ strMsg $ "deleteFeature: Child not present for parent ID " ++ show (f, pid)+                                      else return $ M.adjust (List.delete f) pid m0++parents :: FeatureHier -> Feature -> [Feature]+parents fh f = map parent $ parentIds f+    where parent pid = M.findWithDefault noParent pid $ idToFeature fh+              where noParent = error $ "featParents: Unknown parent id " ++ show pid ++ " for " ++ show f++children :: FeatureHier -> Feature -> [Feature]+children fh f = concatMap chs $ ids f+    where chs fid = M.findWithDefault [] fid $ idToChildren fh++parentsM :: (MonadReader FeatureHier m) => Feature -> m [Feature]+parentsM = asks . flip parents++childrenM :: (MonadReader FeatureHier m) => Feature -> m [Feature]+childrenM = asks . flip children++checkInvariants :: FeatureHier -> [String]+checkInvariants hier+    = concat [ checkAllFeatureIDs, checkAllFeatureParents, checkFeaturesInAll, checkChildrenInAll ]+    where checkAllFeatureIDs = concatMap checkFeatureIDs $ S.toList $ features hier+          checkFeatureIDs f = concatMap checkFeatureID $ ids f+              where checkFeatureID fid = case M.lookup fid $ idToFeature hier of+                                           Nothing -> ["Feature ID " ++ show fid ++ " missing for " ++ show f]+                                           Just f' | f' /= f -> ["Wrong feature " ++ show f' ++ " /= " ++ show f ++ " for ID " ++ show fid]+                                                   | otherwise -> []+          checkAllFeatureParents = []+          checkFeaturesInAll = []+          checkChildrenInAll = []

@@ -0,0 +1,87 @@+{-# OPTIONS_GHC -XFlexibleContexts #-}++module Bio.GFF3.FeatureHierSequences ( FeatureHierSequences, features, sequences+                                     , fromLists, parse+                                     , lookupId, parents, children+                                     , seqData, getSequence, featureSequence+                                     , runGFF, runGFFIO, asksGFF+                                     )+    where ++import Control.Monad.Error+import Control.Monad.Reader+import qualified Data.ByteString.Lazy.Char8 as LBS+import qualified Data.Map as M+import Data.Maybe+import qualified Data.Set as S++import Bio.Sequence.SeqData+import Bio.Sequence.Fasta++import qualified Bio.GFF3.Feature as F+import qualified Bio.GFF3.FeatureHier as FH+import qualified Bio.Location.Location as Loc+import qualified Bio.Location.SeqLocation as SeqLoc+import Bio.Location.OnSeq++data FeatureHierSequences = FeatureHierSequences { featureHier :: !FH.FeatureHier+                                                 , sequenceMap :: !(M.Map SeqName SeqData)+                                                 } deriving (Show)++parse :: (Error e, MonadError e m) => LBS.ByteString -> m FeatureHierSequences+parse str = do (feats, seqlines) <- F.parseWithFasta str+               fromLists feats $ mkSeqs seqlines++fromLists :: (Error e, MonadError e m) => [F.Feature] -> [Sequence] -> m FeatureHierSequences+fromLists feats seqs = let !seqmap = M.fromList $ map seqAssoc seqs+                           !featSeqids = S.fromList $ map F.seqid feats+                           !missingSeqids = featSeqids `S.difference` M.keysSet seqmap+                       in if S.null missingSeqids+                          then liftM (\fh -> FeatureHierSequences fh seqmap) $ FH.fromList feats+                          else throwError $ strMsg $ "Missing sequences for seqids " ++ show (map LBS.unpack $ S.toList missingSeqids)+    where seqAssoc (Seq seqname sequ _) = (seqname, sequ)++features :: FeatureHierSequences -> S.Set F.Feature+features = FH.features . featureHier++sequences :: FeatureHierSequences -> [Sequence]+sequences = M.foldWithKey mkSeqList [] . sequenceMap+    where mkSeqList seqname sequ = ((Seq seqname sequ Nothing) :)++lookupId :: (Error e, MonadError e m) => FeatureHierSequences -> SeqName -> m F.Feature+lookupId fhs = FH.lookupId $ featureHier fhs++children :: FeatureHierSequences -> F.Feature -> [F.Feature]+children = FH.children . featureHier++parents :: FeatureHierSequences -> F.Feature -> [F.Feature]+parents = FH.parents . featureHier++seqData :: (Error e, MonadError e m) => FeatureHierSequences -> SeqLoc.SeqLoc -> m SeqData+seqData fhs (OnSeq chrname loc) = do chrseq <- getSequence fhs chrname+                                     return $ Loc.seqDataPadded chrseq loc++getSequence :: (Error e, MonadError e m) => FeatureHierSequences -> SeqName -> m SeqData+getSequence fhs seqid = maybe seqidNotFound return $ M.lookup seqid $ sequenceMap fhs+    where seqidNotFound = throwError $ strMsg $ "SeqID " ++ show (LBS.unpack seqid) ++ " not found"++featureSequence :: (Error e, MonadError e m) => FeatureHierSequences -> F.Feature -> m Sequence+featureSequence fhs f = do sequ <- seqData fhs $ F.seqLoc f+                           let seqname = fromMaybe (LBS.pack $ SeqLoc.display $ F.seqLoc f) $ listToMaybe $ F.ids f+                           return $ Seq seqname sequ Nothing++catchIOErrors :: IO a -> ErrorT String IO a+catchIOErrors m = ErrorT { runErrorT = liftM Right m `catch` (return . Left . show) }++runGFF :: FilePath -> (ErrorT String (Reader FeatureHierSequences) a) -> ErrorT String IO a+runGFF gffname m = do gff <- catchIOErrors $ LBS.readFile gffname+                      fhs <- parse gff+                      either throwError return $ runReader (runErrorT m) fhs++runGFFIO :: FilePath -> (ErrorT String (ReaderT FeatureHierSequences IO) a) -> ErrorT String IO a+runGFFIO gffname m =  do gff <- catchIOErrors $ LBS.readFile gffname+                         fhs <- parse gff+                         mapErrorT (flip runReaderT fhs) m++asksGFF :: (Error e, MonadError e m, MonadReader FeatureHierSequences m) => (FeatureHierSequences -> a -> m b) -> a -> m b+asksGFF f x = do fhs <- ask; f fhs x

@@ -0,0 +1,94 @@+module Bio.GFF3.SGD ( chromosomes, genes, rRNAs+                    , sortExons+                    , geneSequence, geneSeqLoc, geneCDSes+                    , noncodingSequence, noncodingSeqLoc, noncodingExons+                    , namedSLM, geneCDS_SLM+                    )+    where++import Control.Arrow ((&&&))+import Control.Monad.Error+import Control.Monad.Reader+import qualified Data.ByteString.Lazy.Char8 as LBS+import Data.List+import Data.Maybe+import Data.Ord+import qualified Data.Set as S++import Bio.Sequence.SeqData++import qualified Bio.GFF3.Feature as F+import Bio.GFF3.FeatureHierSequences+-- import qualified Bio.Location.ContigLocation as CLoc+import qualified Bio.Location.Location as Loc+import Bio.Location.OnSeq+import qualified Bio.Location.SeqLocation as SeqLoc+import Bio.Location.Strand+import qualified Bio.Location.SeqLocMap as SLM++-- SGD feature types+cdsType, chromosomeType, geneType, noncodingExonType :: LBS.ByteString+chromosomeType = LBS.pack "chromosome"+cdsType = LBS.pack "CDS"+geneType = LBS.pack "gene"+noncodingExonType = LBS.pack "noncoding_exon"++chromosomes :: FeatureHierSequences -> [F.Feature]+chromosomes = filter isChromosome . S.toList . features+    where isChromosome = (== chromosomeType) . F.ftype++genes :: FeatureHierSequences -> [F.Feature]+genes = filter isGene . S.toList . features+    where isGene = (== geneType) . F.ftype++rRNAs :: FeatureHierSequences -> [F.Feature]+rRNAs = filter isRRNA . S.toList . features+    where isRRNA = any (flip S.member $ S.fromList rRNAFeatureIDs) . F.ids+          rRNAFeatureIDs = map LBS.pack ["RDN25-1", "RDN58-1", "RDN18-1", "RDN5-1"]++geneSequence :: (Error e, MonadError e m) => FeatureHierSequences -> F.Feature -> m Sequence+geneSequence fhs g = do sequ <- geneSeqLoc fhs g >>= SeqLoc.seqData (getSequence fhs)+                        seqname <- maybe (throwError $ strMsg $ "No gene ID for " ++ show g) return $ listToMaybe $ F.ids g+                        return $ Seq seqname sequ Nothing++geneSeqLoc :: (Error e, MonadError e m) => FeatureHierSequences -> F.Feature -> m SeqLoc.SeqLoc+geneSeqLoc fhs g = exonSeqLoc $ geneCDSes fhs g++geneCDSes :: FeatureHierSequences -> F.Feature -> [F.Feature]+geneCDSes fhs g = filter isCDS $ children fhs g+    where isCDS = (== cdsType) . F.ftype          ++noncodingSequence :: (Error e, MonadError e m) => FeatureHierSequences -> F.Feature -> m Sequence+noncodingSequence fhs nc = do sequ <- noncodingSeqLoc fhs nc >>= SeqLoc.seqData (getSequence fhs)+                              seqname <- maybe (throwError $ strMsg $ "No gene ID for " ++ show nc) return $ listToMaybe $ F.ids nc+                              return $ Seq seqname sequ Nothing++noncodingSeqLoc :: (Error e, MonadError e m) => FeatureHierSequences -> F.Feature -> m SeqLoc.SeqLoc+noncodingSeqLoc fhs nc = exonSeqLoc $ noncodingExons fhs nc++noncodingExons :: FeatureHierSequences -> F.Feature -> [F.Feature]+noncodingExons fhs nc = filter isNoncodingExon $ children fhs nc+    where isNoncodingExon = (== noncodingExonType) . F.ftype++sortExons :: (Error e, MonadError e m) => [F.Feature] -> m [F.Feature]+sortExons exons | all fwdStrand exons = return $ sortBy (comparing F.start) exons+                | all rcStrand exons = return $ sortBy (comparing $ negate . F.start) exons+                | otherwise = throwError $ strMsg "sortExons: discordant strands"+    where fwdStrand = (== (Just Fwd)) . F.strand+          rcStrand = (== (Just RevCompl)) . F.strand ++exonSeqLoc :: (Error e, MonadError e m) => [F.Feature] -> m SeqLoc.SeqLoc+exonSeqLoc [] = throwError $ strMsg "exonSeqLoc: empty list of exons"+exonSeqLoc exons@(exon0:exonrest)+    = let seqid = F.seqid exon0+      in if all ((== seqid) . F.seqid) exonrest+         then liftM (OnSeq seqid . Loc.Loc . map F.contigLoc) $ sortExons exons+         else throwError $ strMsg "exonSeqLoc: discordant seqid amongst exons"++namedSLM :: FeatureHierSequences -> SLM.SeqLocMap F.Feature+namedSLM = SLM.fromList . map (F.seqLoc &&& id) . filter isNamed . S.toList . features+    where isNamed = not . null . F.ids++geneCDS_SLM :: (Error e, MonadError e m) => FeatureHierSequences -> m (SLM.SeqLocMap F.Feature)+geneCDS_SLM fhs = liftM SLM.fromList $ mapM cdsLocAndGene $ genes fhs+    where cdsLocAndGene g = liftM (flip (,) g) $ geneSeqLoc fhs g

@@ -0,0 +1,135 @@+-- |Data types for working with locations in a sequence.  Zero-based+--  'Int64' indices are used throughout, to facilitate direct use of+--  indexing functions on 'SeqData'.++module Bio.Location.ContigLocation ( ContigLoc(..), fromStartEnd, fromPosLen+                                   , bounds, startPos, endPos+                                   , slide, extend, posInto, posOutof+                                   , seqData, seqDataPadded, isWithin, overlaps+                                   , display+                                   )+    where ++import Prelude hiding (length)++import Control.Monad.Error+import qualified Data.ByteString.Lazy.Char8 as LBS++import Bio.Sequence.SeqData++import qualified Bio.Location.Position as Pos+import Bio.Location.Strand++-- | Contiguous set of positions in a sequence+data ContigLoc = ContigLoc { offset5 :: !Offset   -- ^ 5' end of region on target sequence, 0-based index+                           , length :: !Offset -- ^ length of region on target sequence+                           , strand :: !Strand -- ^ strand of region+                           } deriving (Eq, Ord, Show)++instance Stranded ContigLoc where+    revCompl (ContigLoc seq5 len str) = ContigLoc seq5 len $ revCompl str++-- | Create a 'ContigLoc' from 0-based starting and ending positions.+--   When 'start' is less than 'end' the position will be on the 'Fwd'+--   'Strand', otherwise it will be on the 'RevCompl' strand.+fromStartEnd :: Offset -> Offset -> ContigLoc+fromStartEnd start end+    | start < end = ContigLoc start (1 + end - start) Fwd+    | otherwise   = ContigLoc end   (1 + start - end) RevCompl++-- | Create a 'ContigLoc' from a Pos.Pos defining the start+-- ('ContigLoc' 5 prime end) position on the sequence and the length.+fromPosLen :: Pos.Pos -> Offset -> ContigLoc+fromPosLen (Pos.Pos off5 Fwd)      len = ContigLoc off5               len Fwd+fromPosLen (Pos.Pos off3 RevCompl) len = ContigLoc (off3 - (len - 1)) len RevCompl++-- | The bounds of a 'ContigLoc', a pair of the lowest and highest+--   sequence indices covered by the region, which ignores the strand+--   of the 'ContigLoc'.  The first element of the pair will always be+--   lower than the second.+bounds :: ContigLoc -> (Offset, Offset)+bounds (ContigLoc seq5 len _) = (seq5, seq5 + len - 1)++-- | 0-based starting (5' in the region orientation) position+startPos :: ContigLoc -> Pos.Pos+startPos (ContigLoc seq5 len str) +    = case str of+        Fwd      -> Pos.Pos seq5             str+        RevCompl -> Pos.Pos (seq5 + len - 1) str++-- | 0-based ending (3' in the region orientation) position+endPos :: ContigLoc -> Pos.Pos+endPos (ContigLoc seq5 len str) +    = case str of+        Fwd      -> Pos.Pos (seq5 + len - 1) str+        RevCompl -> Pos.Pos seq5             str++-- | Move a 'ContigLoc' region by a specified offset+slide :: Offset -> ContigLoc -> ContigLoc+slide dpos (ContigLoc seq5 len str) = ContigLoc (seq5 + dpos) len str++-- | Subsequence 'SeqData' for a 'ContigLoc', provided that the region+--   is entirely within the sequence.+seqData :: (Error e, MonadError e m) => SeqData -> ContigLoc -> m SeqData+seqData sequ (ContigLoc seq5 len str)+    | seq5 < 0 = outOfBounds+    | otherwise = case LBS.take len $ LBS.drop seq5 sequ of+                    fwdseq | LBS.length fwdseq == len -> return $ stranded str fwdseq+                           | otherwise -> outOfBounds+    where outOfBounds = throwError $ strMsg $ "contig seq loc " ++ show (seq5, seq5 + len - 1) ++ " out of SeqData bounds"++-- | Subsequence 'SeqData' for a 'ContigLoc', padded as needed with Ns+seqDataPadded :: SeqData -> ContigLoc -> SeqData+seqDataPadded sequ (ContigLoc seq5 len str) = stranded str fwdseq+    where fwdseq+              | seq5 + len <= 0         = LBS.replicate len 'N'+              | seq5 >= LBS.length sequ = LBS.replicate len 'N'+              | seq5 < 0  = LBS.replicate (negate seq5) 'N' `LBS.append` takePadded (len + seq5) sequ+              | otherwise = takePadded len $ LBS.drop seq5 sequ+          takePadded sublen subsequ+              | sublen <= LBS.length subsequ = LBS.take sublen subsequ+              | otherwise = subsequ `LBS.append` LBS.replicate (sublen - LBS.length subsequ) 'N'++-- | For a 'Pos' and a 'ContigLoc' on the same sequence, find the+--   corresponding 'Pos' relative to the 'ContigLoc', provided it is+--   within the 'ContigLoc'.+posInto :: Pos.Pos -> ContigLoc -> Maybe Pos.Pos+posInto (Pos.Pos pos pStrand) (ContigLoc seq5 len cStrand)+    | pos < seq5 || pos >= seq5 + len = Nothing+    | otherwise = Just $ case cStrand of+                           Fwd      -> Pos.Pos (pos - seq5)              pStrand+                           RevCompl -> Pos.Pos (seq5 + len  - (pos + 1)) (revCompl pStrand)++-- | For a 'Pos' specified relative to a 'ContigLoc', find the+--   corresponding 'Pos' relative to the outer sequence, provided that+--   the 'Pos' is within the bounds of the 'ContigLoc'.+posOutof :: Pos.Pos -> ContigLoc -> Maybe Pos.Pos+posOutof (Pos.Pos pos pStrand) (ContigLoc seq5 len cStrand)+    | pos < 0 || pos >= len = Nothing+    | otherwise = Just $ case cStrand of+                           Fwd      -> Pos.Pos (pos + seq5)             pStrand+                           RevCompl -> Pos.Pos (seq5 + len - (pos + 1)) (revCompl pStrand)++-- | 'ContigLoc' extended on the 5' and 3' ends.+extend :: (Offset, Offset) -> ContigLoc -> ContigLoc+extend (ext5, ext3) (ContigLoc seq5 len str)+    = case str of+        Fwd -> ContigLoc (seq5 - ext5) (len + ext5 + ext3) str+        RevCompl -> ContigLoc (seq5 - ext3) (len + ext5 + ext3) str++-- | For a 'Pos' and a 'ContigLoc' on the same sequence, is the 'Pos'+--   within the 'ContigLoc'.+isWithin :: Pos.Pos -> ContigLoc -> Bool+isWithin (Pos.Pos pos pStrand) (ContigLoc seq5 len cStrand)+    = (pos >= seq5) && (pos < seq5 + len) && (cStrand == pStrand)++-- | For a pair of 'ContigLoc' regions on the same sequence, indicates+--   if they overlap at all.+overlaps :: ContigLoc -> ContigLoc -> Bool+overlaps contig1 contig2+    = case (bounds contig1, bounds contig2) of+        ((low1, high1),(low2, high2)) -> (strand contig1 == strand contig2)+                                         && (low1 <= high2) && (low2 <= high1)++display :: ContigLoc -> String+display cloc = show (Pos.offset $ startPos cloc) ++ "to" ++ show (Pos.offset $ endPos cloc)

@@ -0,0 +1,115 @@+-- | A structure to allow fast lookup of objects whose sequence+--   location lines up with a give position.+module Bio.Location.LocMap (LocMap, mkLocMap, defaultZonesize, fromList+                           , lookupWithin, lookupOverlaps+                           , delete, deleteBy, insert+                           , checkInvariants+                           )+    where ++import qualified Data.IntMap as IM+import qualified Data.IntSet as IS+import Data.List hiding (insert, delete, deleteBy)+import Data.Maybe+import Data.Monoid++import qualified Bio.Location.Location as Loc+import qualified Bio.Location.Position as Pos+import Bio.Sequence.SeqData (Offset)++defaultZonesize :: Offset+defaultZonesize = 65536++-- | Collection mapping a collection of 'Loc' locations, possibly+--   overlapping, binned for efficient lookup by position.+data LocMap a = LocMap !Offset !(IM.IntMap (Loc.Loc, a)) !(IM.IntMap IS.IntSet)++instance Monoid (LocMap a) where+    mempty = mkLocMap defaultZonesize+    mappend lm0 (LocMap _ keyToLoc1 _) = foldl' (\lm (l,x) -> insert l x lm) lm0 $ IM.elems keyToLoc1++-- | Create an empty 'LocMap' with a specified position bin size+mkLocMap :: Offset -> LocMap a+mkLocMap zonesize = LocMap zonesize IM.empty IM.empty++-- | Create a 'LocMap' from an associated list.+fromList :: Offset -> [(Loc.Loc, a)] -> LocMap a+fromList zonesize = foldl' (\lm0 (l,x) -> insert l x lm0) (mkLocMap zonesize)++-- | Add an object with an associated 'Loc' sequence region+insert :: Loc.Loc -> a -> LocMap a -> LocMap a+insert l x (LocMap zonesize keyToLoc zoneKeys) +    = let !key = case IM.maxViewWithKey keyToLoc of+                   Just ((k, _), _) -> k + 1+                   Nothing -> 1+          !newKeyToLoc = IM.insertWithKey duplicateError key (l, x) keyToLoc+          !newZoneKeys = foldl' (insertIntoZone key) zoneKeys $ locZones zonesize l+      in LocMap zonesize newKeyToLoc newZoneKeys+    where insertIntoZone key currZoneKeys z = case IM.lookup z currZoneKeys of+                                                Nothing -> IM.insert z (IS.singleton key) currZoneKeys+                                                Just currZoneKeySet -> IM.insert z (IS.insert key currZoneKeySet) currZoneKeys+          duplicateError k (l1, _) (l2, _) = error $ "LocMap: key " ++ show k ++ " duplicated: " ++ show (l1, l2)++-- | Find the (possibly empty) list of sequence regions and associated+--   objects that contain a 'Pos' position, in the sense of 'withinLoc'+lookupWithin :: Pos.Pos -> LocMap a -> [(Loc.Loc, a)]+lookupWithin pos (LocMap zonesize keyToLoc zoneKeys) +    = let !zoneKeySet = IM.findWithDefault IS.empty (posZone zonesize pos) zoneKeys+      in filter ((Loc.isWithin pos) . fst) $ keySetLocs keyToLoc zoneKeySet++-- | Find the (possibly empty) list of sequence regions and associated+--   objects that overlap a 'Loc' region, in the sense of 'overlapsLoc'+lookupOverlaps :: Loc.Loc -> LocMap a -> [(Loc.Loc, a)]+lookupOverlaps loc (LocMap zonesize keyToLoc zoneKeys)+    = let !zonesKeySet = IS.unions $ map (\z -> IM.findWithDefault IS.empty z zoneKeys) $ locZones zonesize loc+      in filter ((Loc.overlaps loc) . fst) $ keySetLocs keyToLoc zonesKeySet++-- | Remove a region / object association from the map, if it is+-- present.  If it is present multiple times, only the first+-- occurrence will be deleted.+delete :: (Eq a) => (Loc.Loc, a) -> LocMap a -> LocMap a+delete target = deleteBy (== target)++-- | Remove the first region / object association satisfying a+-- predicate function.+deleteBy :: ((Loc.Loc, a) -> Bool) -> LocMap a -> LocMap a+deleteBy isTarget lm0@(LocMap zonesize keyToLoc zoneKeys) +    = case filter (isTarget . snd) $ IM.toList keyToLoc of+        [] -> lm0+        ((key,(l,_)):_) -> let !newKeyToLoc = IM.delete key keyToLoc+                               !newZoneKeys = foldl' (deleteFromZone key) zoneKeys $ locZones zonesize l+                           in LocMap zonesize newKeyToLoc newZoneKeys+    where deleteFromZone key currZoneKeys z = let !currZoneKeySet = IM.findWithDefault missingZone z currZoneKeys +                                              in IM.insert z (IS.delete key currZoneKeySet) currZoneKeys+              where missingZone = error $ "LocMap deleteBy: empty keyset for zone " ++ show z++posZone :: Offset -> Pos.Pos -> Int+posZone zonesize pos = fromIntegral $ (Pos.offset pos) `div` zonesize++locZones :: Offset -> Loc.Loc -> [Int]+locZones zonesize loc = let !(pmin, pmax) = Loc.bounds loc+                            !zmin = fromIntegral $ pmin `div` zonesize+                            !zmax = fromIntegral $ pmax `div` zonesize+                        in [zmin..zmax]++keySetLocs :: (IM.IntMap (Loc.Loc, a)) -> IS.IntSet -> [(Loc.Loc, a)]+keySetLocs keyToLoc = map keyLoc . IS.toList+    where keyLoc key = IM.findWithDefault unknownKey key keyToLoc+              where unknownKey = error $ "LocMap: unknown key " ++ show key++checkInvariants :: LocMap a -> [String]+checkInvariants (LocMap zonesize keyToLoc zoneKeys)+    = concat $ checkAllKeys : map checkKeyLoc (IM.toList keyToLoc)+    where checkAllKeys = let !allZoneKeys = IS.unions $ IM.elems zoneKeys +                             !allLocKeys = IM.keysSet keyToLoc+                             !missingKeys = IS.toAscList $ allZoneKeys `IS.difference` allLocKeys+                             !extraKeys = IS.toAscList $ allLocKeys `IS.difference` allZoneKeys+                         in concat [ map (\key -> "Missing key " ++ show key) missingKeys+                                   , map (\key -> "Extra key " ++ show key) extraKeys]+          checkKeyLoc (key, (loc, _)) = let !keyZoneSet = IM.keysSet $ IM.filter (\keyset -> key `IS.member` keyset) zoneKeys+                                            !locZoneSet = IS.fromList $ locZones zonesize loc+                                            !missingZones = IS.toAscList $ locZoneSet `IS.difference` keyZoneSet+                                            !extraZones = IS.toAscList $ keyZoneSet `IS.difference` locZoneSet+                                        in concat [ map (\zone -> "Missing zone " ++ show zone ++ " for " ++ show (key, loc)) missingZones+                                                  , map (\zone -> "Extra zone " ++ show zone ++ " for " ++ show (key, loc)) extraZones]+

@@ -0,0 +1,119 @@+-- |Data types for working with locations in a sequence.  Zero-based+--  'Int64' indices are used throughout, to facilitate direct use of+--  indexing functions on 'SeqData'.++module Bio.Location.Location ( Loc(..), bounds, length, startPos, endPos+                             , extend, posInto, posOutof, isWithin, overlaps+                             , seqData, seqDataPadded+                             , display+                             )+    where ++import Prelude hiding (length)++import Control.Arrow ((***))+import Control.Monad.Error+import qualified Data.ByteString.Lazy.Char8 as LBS+import Data.List (intercalate)+--import Data.Maybe++import qualified Bio.Location.ContigLocation as CLoc+import qualified Bio.Location.Position as Pos+import Bio.Location.Strand+import Bio.Sequence.SeqData++-- | General (disjoint) sequence region consisting of a concatenated+--   set of contiguous regions+newtype Loc = Loc [CLoc.ContigLoc] deriving (Eq, Ord, Show)++instance Stranded Loc where+    revCompl (Loc contigs) = Loc $ reverse $ map revCompl contigs++-- | Length of the region+length :: Loc -> Offset+length (Loc contigs) = sum $ map CLoc.length contigs++-- | The bounds of a 'Loc', consisting of the lowest & highest+--   sequence indices lying within the region.  The first element of+--   the pair will always be lower than the second.+bounds :: Loc -> (Offset, Offset)+bounds (Loc []) = error "locBounds on zero-contig Loc"+bounds (Loc contigs) = (minimum *** maximum) $ unzip $ map CLoc.bounds contigs++-- | 0-based starting (5' in the region orientation) offset of the+--   region on its sequence.+startPos :: Loc -> Pos.Pos+startPos (Loc [])      = error "startPos: zero-contig Loc"+startPos (Loc contigs) = CLoc.startPos $ head contigs++-- | 0-based ending (3' in the region orientation) offset of the+--   region on its sequence.+endPos :: Loc -> Pos.Pos+endPos (Loc [])      = error "endPos: zero-contig Loc"+endPos (Loc contigs) = CLoc.endPos $ last contigs++-- | Subsequence 'SeqData' for a 'Loc', provided that the region is+--   entirely within the sequence.+seqData :: (Error e, MonadError e m) => SeqData -> Loc -> m SeqData+seqData sequ (Loc contigs)+    = liftM LBS.concat $ mapM (CLoc.seqData sequ) contigs++seqDataPadded :: SeqData -> Loc -> SeqData+seqDataPadded sequ (Loc contigs)+    = LBS.concat $ map (CLoc.seqDataPadded sequ) contigs++-- | For a 'Pos' and a 'Loc' region on the same sequence, find the+--   corresponding 'Pos' relative to the region, if the 'Pos' is+--   within the region.  If the 'Loc' region has redundant positions+--   for a given sequence position, the first is returned.+posInto :: Pos.Pos -> Loc -> Maybe Pos.Pos+posInto seqpos (Loc contigs) = posIntoContigs seqpos contigs++posIntoContigs :: Pos.Pos -> [CLoc.ContigLoc] -> Maybe Pos.Pos+posIntoContigs _ [] = Nothing+posIntoContigs seqpos (contig@(CLoc.ContigLoc _ len _):rest)+    = case CLoc.posInto seqpos contig of+        just@(Just _) -> just+        Nothing -> liftM (flip Pos.slide len) $ posIntoContigs seqpos rest++-- | For a 'Loc' region on a sequence and a 'Pos' relative to the+--   region, find the corresponding 'Pos' on the sequence, provided+--   that the position is within the bounds of the region.+posOutof :: Pos.Pos -> Loc -> Maybe Pos.Pos+posOutof pos (Loc contigs) = posOutofContigs pos contigs++posOutofContigs :: Pos.Pos -> [CLoc.ContigLoc] -> Maybe Pos.Pos+posOutofContigs _ [] = Nothing+posOutofContigs seqpos (contig@(CLoc.ContigLoc _ len _):rest)+    = case CLoc.posOutof seqpos contig of+        just@(Just _) -> just+        Nothing -> posOutofContigs (Pos.slide seqpos $ negate len) rest++-- | Extend a 'Loc' region by incorporating contigous nucleotide+--   regions of the specified lengths on the 5' and 3' ends+extend :: (Offset, Offset) -- ^ (5' extension, 3' extension)+       -> Loc -> Loc+extend _ (Loc []) = error "extendLoc on zero-contig Loc"+extend (ext5, ext3) (Loc contigs) = Loc $ case extendContigs3 contigs of+                                            [] -> error "Empty contigs after extendContigs3"+                                            (cfirst:crest) -> (CLoc.extend (ext5, 0) cfirst):crest+    where extendContigs3 [] = error "Empty contigs in extendContigs3"+          extendContigs3 [clast] = [CLoc.extend (0, ext3) clast]+          extendContigs3 (contig:crest) = contig : extendContigs3 crest++-- | For a 'Pos' and a 'Loc' on the same sequence, does the position+--   fall within the 'Loc' region?+isWithin :: Pos.Pos -> Loc -> Bool+isWithin seqpos (Loc contigs) = or $ map (CLoc.isWithin seqpos) contigs++overlappingContigs :: Loc -> Loc -> [(CLoc.ContigLoc, CLoc.ContigLoc)]+overlappingContigs (Loc contigs1) (Loc contigs2) +    = filter (uncurry CLoc.overlaps) [(c1, c2) | c1 <- contigs1, c2 <- contigs2 ]++-- | For a pair of 'Loc' regions on the same sequence, do they overlap+--   at all?+overlaps :: Loc -> Loc -> Bool+overlaps l1 l2 = not $ null $ overlappingContigs l1 l2++display :: Loc -> String+display (Loc contigs) = intercalate ";" $ map CLoc.display contigs

@@ -0,0 +1,48 @@+module Bio.Location.OnSeq ( SeqName +                          , OnSeq(..), withSeqData, andSameSeq, onSameSeq+                          , OnSeqs, perSeq, perSeqUpdate, withNameAndSeq+                          )+    where ++import Control.Monad.Error+import qualified Data.ByteString.Lazy.Char8 as LBS+import Data.List+import qualified Data.Map as M+import Data.Monoid++import Bio.Sequence.SeqData++type SeqName = SeqData++data OnSeq a = OnSeq { onSeqName :: !SeqName+                     , onSeqObj :: !a+                     } deriving (Eq, Ord, Show)++instance Functor OnSeq where+    fmap f (OnSeq seqname x) = OnSeq seqname (f x)++withSeqData :: (Error e, MonadError e m) => (SeqData -> a -> m b) -> (SeqName -> m SeqData) -> OnSeq a -> m b+withSeqData f lookupSeq (OnSeq seqname x) = lookupSeq seqname >>= flip f x++andSameSeq :: (a -> b -> Bool) -> OnSeq a -> OnSeq b -> Bool+andSameSeq f (OnSeq xname x) (OnSeq yname y) | xname == yname = f x y+                                             | otherwise = False++onSameSeq :: (Monad m) => (a -> b -> m c) -> OnSeq a -> OnSeq b -> m c+onSameSeq f (OnSeq xname x) (OnSeq yname y) +    | xname == yname = f x y+    | otherwise      = fail $ "onSameSeq: " ++ show (LBS.unpack xname) ++ " /= " ++ show (LBS.unpack yname)+++type OnSeqs a = M.Map SeqName a++perSeq :: (Monoid b) => (a -> b -> c) -> OnSeq a -> OnSeqs b -> c+perSeq f (OnSeq seqname x) = f x . M.findWithDefault mempty seqname++perSeqUpdate :: (Monoid b) => (a -> b -> b) -> OnSeq a -> OnSeqs b -> OnSeqs b+perSeqUpdate upd onseq@(OnSeq seqname _) seqmap0 = M.insert seqname (perSeq upd onseq seqmap0) seqmap0++withNameAndSeq :: (Monad m) => (SeqName -> a -> b -> m c) -> OnSeq a -> OnSeqs b -> m c+withNameAndSeq f (OnSeq seqname x) = mylookup seqname >=> f seqname x+    where mylookup k = maybe nameNotFound return . M.lookup k+              where nameNotFound = fail $ "withNameAndSeq: sequence " ++ show (LBS.unpack k) ++ " not found"

@@ -0,0 +1,41 @@+-- | Positions on a sequence.  Zero-based 'Int64' indices are used+-- throughout, to facilitate direct use of indexing functions on+-- 'SeqData'.++module Bio.Location.Position ( Pos(..), slide, seqNt, seqNtPadded, display )+    where ++import Control.Monad.Error+import qualified Data.ByteString.Lazy.Char8 as LBS+import Data.Ix (Ix)+import Data.List++import Bio.Location.Strand+import Bio.Sequence.SeqData++-- | Position in a sequence+data Pos = Pos { offset :: !Offset -- ^ 0-based index of the position+               , strand :: !Strand -- ^ Optional strand of the position+               }+              deriving (Eq, Ord, Show, Read, Ix)++instance Stranded Pos where+    revCompl (Pos off str) = Pos off $ revCompl str++-- | Slide a position by an offset+slide :: Pos -> Offset -> Pos+slide (Pos off str) doff = Pos (off + doff) str++seqNt :: (Error e, MonadError e m) => SeqData -> Pos -> m Char+seqNt sequ (Pos off str) | off >= 0 && off < LBS.length sequ = return $ stranded str $ sequ `LBS.index` off+                         | otherwise = throwError $ strMsg $ "position " ++ show off ++ " out of SeqData bounds"++seqNtPadded :: SeqData -> Pos -> Char+seqNtPadded sequ (Pos off str) | off >= 0 && off < LBS.length sequ = stranded str $ sequ `LBS.index` off+                               | otherwise = 'N'++display :: Pos -> String+display (Pos off str) = show off ++ displayStrand+    where displayStrand = case str of+                            Fwd -> "(+)"+                            RevCompl -> "(-)"

@@ -0,0 +1,33 @@+module Bio.Location.SeqLocMap ( SeqLocMap+                              , empty+                              , fromList, insert+                              , lookupWithin, lookupOverlaps+                              )+    where ++import Control.Arrow (first)+import Data.List hiding (insert)+import qualified Data.Map as M++import Bio.Location.OnSeq+import qualified Bio.Location.LocMap as LM+import qualified Bio.Location.SeqLocation as SeqLoc++type SeqLocMap a = OnSeqs (LM.LocMap a)++empty :: SeqLocMap a+empty = M.empty++fromList :: [(SeqLoc.SeqLoc, a)] -> SeqLocMap a+fromList = foldl' (\lm (sl,x) -> insert sl x lm) M.empty++insert :: SeqLoc.SeqLoc -> a -> SeqLocMap a -> SeqLocMap a+insert sloc x = perSeqUpdate (\loc locmap -> LM.insert loc x locmap) sloc++lookupWithin :: SeqLoc.SeqPos -> SeqLocMap a -> [(SeqLoc.SeqLoc, a)]+lookupWithin = withNameAndSeq namedLookupWithin+    where namedLookupWithin seqname pos = map (first $ OnSeq seqname) . LM.lookupWithin pos++lookupOverlaps :: SeqLoc.SeqLoc -> SeqLocMap a -> [(SeqLoc.SeqLoc, a)]+lookupOverlaps = withNameAndSeq namedLookupOverlaps+    where namedLookupOverlaps seqname loc = map (first $ OnSeq seqname) . LM.lookupOverlaps loc

@@ -0,0 +1,46 @@+module Bio.Location.SeqLocation ( SeqPos, displaySeqPos+                                , ContigSeqLoc, withinContigSeqLoc, displayContigSeqLoc +                                , SeqLoc+                                , isWithin, overlaps, seqData+                                , display+                                )+    where ++import Control.Monad.Error+import qualified Data.ByteString.Lazy.Char8 as LBS+--import Data.List++import qualified Bio.Location.ContigLocation as CLoc+import qualified Bio.Location.Location as Loc+import Bio.Location.OnSeq+import qualified Bio.Location.Position as Pos+import Bio.Sequence.SeqData++type SeqPos = OnSeq Pos.Pos++displaySeqPos :: SeqPos -> String+displaySeqPos (OnSeq refname pos) = LBS.unpack refname ++ "@" ++ Pos.display pos++type ContigSeqLoc = OnSeq CLoc.ContigLoc++withinContigSeqLoc :: SeqPos -> ContigSeqLoc -> Bool+withinContigSeqLoc = andSameSeq CLoc.isWithin++displayContigSeqLoc :: ContigSeqLoc -> String+displayContigSeqLoc (OnSeq refname cloc) = LBS.unpack refname ++ "@" ++ CLoc.display cloc++type SeqLoc = OnSeq Loc.Loc++isWithin :: SeqPos -> SeqLoc -> Bool+isWithin = andSameSeq Loc.isWithin++overlaps :: SeqLoc -> SeqLoc -> Bool+overlaps = andSameSeq Loc.overlaps++seqData :: (Error e, MonadError e m) => (SeqName -> m SeqData) -> SeqLoc -> m SeqData+seqData = withSeqData Loc.seqData++display :: SeqLoc -> String+display (OnSeq refname loc) = LBS.unpack refname ++ "@" ++ Loc.display loc++

@@ -0,0 +1,33 @@+module Bio.Location.Strand ( Strand(..)+                           , Stranded(..), stranded+                           )+    where ++import qualified Data.ByteString.Lazy.Char8 as LBS+import Data.Ix (Ix)++import Bio.Sequence.SeqData++-- | Sequence strand+data Strand = Fwd | RevCompl deriving (Eq, Ord, Show, Read, Bounded, Enum, Ix)++-- | Anything, such as a location or a sequence, which lies on a+-- strand and can thus be reverse complemented.+class Stranded s where+    revCompl :: s -> s++stranded :: (Stranded s) => Strand -> s -> s+stranded Fwd      = id+stranded RevCompl = revCompl++instance Stranded Strand where+    revCompl Fwd      = RevCompl+    revCompl RevCompl = Fwd++instance Stranded Char where+    revCompl = compl++-- Avoid using -XTypeSynonymInstances+instance Stranded LBS.ByteString where+    revCompl = LBS.reverse . LBS.map compl+

@@ -29,6 +29,11 @@     , readFastaQual     , writeFastaQual, hWriteFastaQual +    -- ** The FastQ format ("Bio.Sequence.FastQ") +    -- Combines sequence data and quality in one file.+    -- Warning: Solexa uses a different formula for the quality values!+    , readFastQ, writeFastQ, hReadFastQ, hWriteFastQ+     -- ** The phd file format ("Bio.Sequence.Phd")     -- | These contain base (nucleotide) calling information,     --   and are generated by @phred@.@@ -53,6 +58,7 @@  -- file formats import Bio.Sequence.Fasta+import Bio.Sequence.FastQ import Bio.Sequence.Phd import Bio.Sequence.TwoBit 

@@ -0,0 +1,70 @@+{- |+   Module: Bio.Sequence.FastQ++   Support the FastQ format that combines sequence and quality. See:++   * <http://www.bioperl.org/wiki/FASTQ_sequence_format>++   Of course, this is yet another vaguely defined pseudo-standard with +   conflicting definitions.  Of course Solexa had to go and invent a different, +   but indistinguishably so, way to do it:++   * <http://www.bcgsc.ca/pipermail/ssrformat/2007-March/000137.html>++   * <http://maq.sourceforge.net/fastq.shtml>++   Currently, we only support the non-Solexa FastQ, adding\/subtracting 33 for +   the quality values.+-}++module Bio.Sequence.FastQ +    (+     -- * Reading FastQ+    readFastQ, hReadFastQ, parse+     -- * Writing FastQ+    , writeFastQ, hWriteFastQ, unparse+    ) where++import System.IO+import qualified Data.ByteString.Lazy.Char8 as B+import qualified Data.ByteString.Lazy as BB+import Data.List (unfoldr)++import Bio.Sequence.SeqData++readFastQ :: FilePath -> IO [Sequence]+readFastQ f = (go . B.lines) `fmap` B.readFile f ++hReadFastQ :: Handle -> IO [Sequence]+hReadFastQ h = (go . B.lines) `fmap` B.hGetContents h++go :: [B.ByteString] -> [Sequence]+go = map (either error id) . unfoldr parse++-- | Parse one FastQ entry, suitable for using in 'unfoldr' over+--   'B.lines' from a file+parse :: [B.ByteString] -> Maybe (Either String Sequence, [B.ByteString])+parse (h1:sd:h2:sq:rest) = +    case (B.uncons h1,B.uncons h2) of+      (Just ('@',h1name), Just ('+',h2name))+          | h1name == h2name -> Just (Right $ Seq h1name sd (Just (BB.map (subtract 33) sq)), rest)+          | otherwise        -> Just (Left $ "Bio.Sequence.FastQ: name mismatch:" ++ showStanza, rest)+      _                      -> Just (Left $ "Bio.Sequence.FastQ: illegal FastQ format:" ++ showStanza, rest)+    where showStanza = unlines $ map B.unpack [ h1, sd, h2, sq ]+parse [] = Nothing+parse fs = let showStanza = unlines (map B.unpack fs)+               err = Left $ "Bio.Sequence.FastQ: illegal number of lines in FastQ format: " ++ showStanza+           in Just (err, [])++writeFastQ :: FilePath -> [Sequence] -> IO ()+writeFastQ f = B.writeFile f . B.concat . map unparse++hWriteFastQ :: Handle -> [Sequence] -> IO ()+hWriteFastQ h = B.hPut h . B.concat . map unparse++-- helper function for writing+unparse :: Sequence -> B.ByteString+unparse (Seq h sd (Just sq)) = +    B.unlines [B.cons '@' h, sd, B.cons '+' h, BB.map (+33) sq]+unparse (Seq h _ Nothing) = error ("Bio.Sequence.FastQ: sequence " ++ show (B.unpack h) +                                   ++" doesn not have quality data!")

@@ -19,6 +19,8 @@     , countSeqs     -- * Helper function for reading your own sequences     , mkSeqs+    -- * Data structures+    , Qual ) where  @@ -27,10 +29,8 @@ import Data.Int import Data.Maybe import System.IO-import System.IO.Unsafe import Control.Monad -import qualified Data.ByteString.Char8 as BS import qualified Data.ByteString.Lazy.Char8 as B import qualified Data.ByteString.Lazy as BB import Data.ByteString.Lazy.Char8 (ByteString)@@ -58,9 +58,7 @@  -- | Lazily read sequences from a FASTA-formatted file readFasta :: FilePath -> IO [Sequence]-readFasta f = do-  ls <- getLines f-  return (mkSeqs ls)+readFasta f = (mkSeqs . B.lines) `fmap` B.readFile f  -- | Write sequences to a FASTA-formatted file. --   Line length is 60.@@ -72,9 +70,7 @@  -- | Read quality data for sequences to a file. readQual :: FilePath -> IO [Sequence]-readQual f = do-  ls <- getLines f-  return (mkQual ls)+readQual f = (mkQual . B.lines) `fmap` B.readFile f  -- | Write quality data for sequences to a file. writeQual :: FilePath -> [Sequence] -> IO ()@@ -89,7 +85,6 @@ readFastaQual  s q = do   ss <- readFasta s   qs <- readQual q-  return ss   -- warning: assumes correct qual file!   let mkseq s1@(Seq x y _) s2@(Seq _ _ (Just z))           | seqlabel s1 /= seqlabel s2 = error ("mismatching sequence and quality: "@@ -98,7 +93,10 @@                                                 ++ 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+      zipWith' f (x:xs) (y:ys) = f x y : zipWith' f xs ys+      zipWith' _ [] [] = []+      zipWith' _ _ _ = error "readFastaQual: Unbalanced sequence and quality"+  return $ zipWith' mkseq ss qs  -- | Write sequence and quality data simulatnously --   This may be more laziness-friendly.@@ -116,9 +114,7 @@  -- | Lazily read sequence from handle hReadFasta :: Handle -> IO [Sequence]-hReadFasta h = do-  ls <- hGetLines h-  return (mkSeqs ls)+hReadFasta h = (mkSeqs . B.lines) `fmap` B.hGetContents h  -- | Write sequences in FASTA format to a handle. hWriteFasta :: Handle -> [Sequence] -> IO ()@@ -147,31 +143,6 @@       qs = B.pack . unwords . map show . BB.unpack   mapM_ ((\l' -> B.hPut h l' >> B.hPut h (B.pack "\n")) . qs) qls wQual _ (Seq _ _ Nothing) = return ()---- ByteString operations--- These aren't (or possible weren't) provided by the FPS library.--- Implement line-based IO (in retrospect it'd be simpler--- and better to just use 'lines'.)---- lazily read lines from file-getLines :: FilePath -> IO [ByteString]-getLines f = do-  h <- openFile f ReadMode-  hGetLines' (hClose h) h---- lazily read lines from handle-hGetLines :: Handle -> IO [ByteString]-hGetLines = hGetLines' (return ())---- add an optional handle-closing parameter-hGetLines' :: IO () -> Handle -> IO [ByteString]-hGetLines' c h = do-  e <- hIsEOF h-  if e then c >> return []-       else do l' <- BS.hGetLine h-               let l = B.fromChunks $ if BS.null l' then [] else  [l']-               ls <- unsafeInterleaveIO $ hGetLines' c h-               return (l:ls)  -- | Convert a list of FASTA-formatted lines into a list of sequences. --   Blank lines are ignored.

@@ -16,3 +16,33 @@  module Bio.Sequence.KEGG where +import Bio.Sequence.GeneOntology (UniProtAcc)++import qualified Data.ByteString.Lazy.Char8 as B+import Data.ByteString.Lazy.Char8 (ByteString)+import Control.Arrow ((&&&))++-- | Most KEGG files that contain associations, have one association per line,+--   consisting of two items separated by whitespace.  This is a generalized reader+--   function.+genReadKegg :: FilePath -> IO [(ByteString,ByteString)]+genReadKegg f = return . map (((!!0) &&& (!!1)) . B.words) . B.lines =<< B.readFile f++newtype KO = KO ByteString+instance Show KO where show (KO x) = B.unpack x++-- | Convert UniProt IDs (up:xxxxxx) to the "UniProtAcc" type.+decodeUP :: ByteString -> UniProtAcc+decodeUP = removePrefix "up" "UniProt accession" id++-- | Convert KO IDs (ko:xxxxx) to the "KO" data type.+decodeKO :: ByteString -> KO+decodeKO = removePrefix "ko" "KEGG KO id" KO++-- | KEGG uses strings with an identifying prefix for IDs. This helper function checks+--   and removes prefix to construct native values.+removePrefix :: String -> String ->  (ByteString -> a) -> ByteString -> a+removePrefix pfx err conv bs+    | (B.pack (pfx++":") `B.isPrefixOf` bs) = conv $ B.drop (fromIntegral (length pfx+1)) bs+    | otherwise                             = error ("Can't parse as "++err ++":\n"+                                                        ++B.unpack bs)

@@ -0,0 +1,247 @@+{- | Read (and write?) the SFF file format used by+   Roche\/454 sequencing to store flowgram data.++   A flowgram is a series of values (intensities) representing homopolymer runs of+   A,G,C, and T in a fixed cycle, and usually displayed as a histogram.++   The Staden Package contains an io_lib, with a C routine for parsing this format.+   According to comments in the sources, the io_lib implementation is based on a file+   called getsff.c, which I've been unable to track down.++   It is believed that all values are stored big endian.+-}++module Bio.Sequence.SFF ( SFF(..), CommonHeader(..)+                        , ReadHeader(..), ReadBlock(..)+                        , readSFF, writeSFF+                        , sffToSequence+                        , test, convert+                        , Flow, Qual, Index+                        ) where++import Bio.Sequence.SeqData++import Data.Int+import qualified Data.ByteString.Lazy as LB+import qualified Data.ByteString as B+import qualified Data.ByteString.Char8 as BC+import Data.ByteString (ByteString)+import Control.Monad (when,replicateM,replicateM_)++import Data.List (intersperse)+import Data.Binary+import Data.Binary.Get (getByteString)+import qualified Data.Binary.Get as G+import Data.Binary.Put (putByteString)++import Text.Printf (printf)++-- | The type of flowgram value+type Flow = Int16+type Index = Word8++-- Global variables holding static information+magic :: Int32+magic = 0x2e736666++versions :: [Int32]+versions = [1]+++readSFF :: FilePath -> IO SFF+readSFF f = return . decode =<< LB.readFile f++sffToSequence :: SFF -> [Sequence]+sffToSequence (SFF ch rs) = map r2s rs+    where r2s r = clip (read_name $ read_header r, bases r, quality r)+          lb x = LB.fromChunks [x]+          clip (n,s,q) = let (k,s2) = B.splitAt (fromIntegral $ key_length ch) s+                         in if k==key ch then Seq (lb n) (lb s2) (Just $ lb $ B.drop (fromIntegral $ key_length ch) q)+                                else error ("Couldn't match key in sequence "++BC.unpack n++" ("++BC.unpack k++" vs. "++BC.unpack (key ch)++")!")++writeSFF :: FilePath -> SFF -> IO ()+writeSFF = encodeFile++-- --------------------------------------------------+-- | test serialization by output'ing the header and first two reads +--   in an SFF, and the same after a decode + encode cycle.+test :: FilePath -> IO ()+test file = do +  (SFF h rs) <- readSFF file +  let sff = (SFF h (take 2 rs))+  putStrLn $ show $ sff+  putStrLn ""+  putStrLn $ show $ (decode $ encode sff :: SFF)++-- --------------------------------------------------+-- | Convert a file by decoding it and re-encoding it+--   This will lose the index (which isn't really necessary)+convert :: FilePath -> IO ()+convert file = writeSFF (file++".out") =<< readSFF file++-- | Generalized function for padding+pad :: Integral a => a -> Put+pad x = replicateM_ (fromIntegral x) (put zero) where zero = 0 :: Word8 ++-- | Generalized function to skip padding+skip :: Integral a => a -> Get ()+skip = G.skip . fromIntegral++-- | The data structure storing the contents of an SFF file (modulo the index)+data SFF = SFF !CommonHeader [ReadBlock]++instance Show SFF where +    show (SFF h rs) = (show h ++ "Read Blocks:\n\n" ++ concatMap show rs)++instance Binary SFF where+    get = do+      -- Parse CommonHeader+      chead <- get+      -- Parse ReadHeaders+      rds <- replicateM (fromIntegral (num_reads chead))+                                   (do +                                      rh <- get :: Get ReadHeader+                                      let nb = fromIntegral $ num_bases rh+                                          fl = fromIntegral $ flow_length chead+                                      fg <- getByteString (2*fl)+                                      fi <- getByteString nb+                                      bs <- getByteString nb+                                      qty <- getByteString nb+                                      let l = (fl*2+nb*3) `mod` 8+                                      when (l > 0) (skip (8-l))+                                      return (ReadBlock rh (decodeArray fg) fi bs qty)+                                   )+      return (SFF chead rds)++    put (SFF hd rds) = do+      put hd+      mapM_ (putRB (fromIntegral $ flow_length hd)) rds++-- | A ReadBlock can't be an instance of Binary directly, since it depends on+--   information from the CommonHeader.+putRB :: Int -> ReadBlock -> Put+putRB fl rb = do+  put (read_header rb)+  mapM_ put (flowgram rb)+  putByteString (flow_index rb)+  putByteString (bases rb)+  putByteString (quality rb)+  let nb = fromIntegral $ num_bases $ read_header rb+      l = (fl*2+nb*3) `mod` 8+  when (l > 0) (pad (8-l))++-- | What the name and type says.+decodeArray :: ByteString -> [Flow]+decodeArray = dec . map fromIntegral . B.unpack +    where dec (d1:d2:rest) = d1*256+d2 : dec rest+          dec [] = []+          dec _  = error "odd flowgram length?!"++-- ----------------------------------------------------------+-- | SFF has a 31-byte common header+--   Todo: remove items that are derivable (counters, magic, etc)+--   cheader_lenght points to the first read header.+--   Also, the format is open to having the index anywhere between reads,+--   we should really keep count and check for each read.  In practice, it+--   seems to be places after the reads.+--   +--   The following two fields are considered part of the header, but as+--   they are static, they are not part of the data structure+--     magic   :: Word32   -- ^ 0x2e736666, i.e. the string ".sff"+--     version :: Word32   -- ^ 0x00000001++data CommonHeader = CommonHeader {+          index_offset                            :: Int64    -- ^ Points to a text(?) section+        , index_length, num_reads                 :: Int32+        , key_length, flow_length                 :: Int16+        , flowgram_fmt                            :: Word8+        , flow, key                               :: ByteString +        }++instance Show CommonHeader where+    show (CommonHeader io il nr kl fl fmt f k) =+        "Common Header:\n\n" ++ (unlines $ map ("    "++) +                                 ["index_off:\t"++show io ++"\tindex_len:\t"++show il+                                 ,"num_reads:\t"++show nr+                                 ,"key_len:\t"  ++show kl ++ "\tflow_len:\t"++show fl+                                 ,"format\t:"   ++show fmt+                                 ,"flow\t:"     ++BC.unpack f+                                 ,"key\t:"      ++BC.unpack k+                                 , ""+                                 ])++instance Binary CommonHeader where+    get = do { m <- get ; when (m /= magic)   $ error (printf "Incorrect magic number - got %8x, expected %8x" m magic)+             ; v <- get ; when (not (v `elem` versions)) $ error (printf "Unexpected version - got %d, supported are: %s" v (unwords $ map show versions))+             ; io <- get ; ixl <- get ; nrd <- get+             ; chl <- get ; kl <- get ; fl <- get ; fmt <- get+             ; fw <- getByteString (fromIntegral fl)+             ; k  <- getByteString (fromIntegral kl)+             ; skip (chl-(31+fl+kl)) -- skip to boundary+             ; return (CommonHeader io ixl nrd kl fl fmt fw k)+             }++    put ch = let CommonHeader io il nr kl fl fmt f k = ch { index_offset = 0 } in+        do { let cl = 31+fl+kl+                 l = cl `mod` 8+                 padding = if l > 0 then 8-l else 0+           ; put magic; put (last versions); put io; put il; put nr; put (cl+padding); put kl; put fl; put fmt+           ; putByteString f; putByteString k+           ; pad padding -- skip to boundary+           }++-- ---------------------------------------------------------- +-- | Each Read has a fixed read header+data ReadHeader = ReadHeader {+      name_length                           :: Int16+    , num_bases                             :: Int32+    , clip_qual_left, clip_qual_right+    , clip_adapter_left, clip_adapter_right :: Int16+    , read_name                             :: ByteString+}++instance Show ReadHeader where+    show (ReadHeader nl nb cql cqr cal car rn) =+        ("    Read Header:\n" ++) $ unlines $ map ("        "++) +                    [ "name_len:\t"++show nl, "num_bases:\t"++show nb+                    , "clip_qual:\t"++show cql++"..."++show cqr+                    , "clip_adap:\t"++show cal++"..."++show car+                    , "read name:\t"++BC.unpack rn+                    , "" +                    ]++instance Binary ReadHeader where+    get = do+      { rhl <- get; nl <- get; nb <- get+      ; cql <- get; cqr <- get ; cal <- get ; car <- get+      ; n <- getByteString (fromIntegral nl)+      ; skip (rhl - (16+ nl))+      ; return (ReadHeader nl nb cql cqr cal car n)+      }+    put (ReadHeader nl nb cql cqr cal car rn) = +        do { let rl = 16+nl+                 l = rl `mod` 8+                 padding = if l > 0 then 8-l else 0+           ; put (rl+padding); put nl; put nb; put cql; put cqr; put cal; put car+           ; putByteString rn +           ; pad padding+           }++-- ----------------------------------------------------------+-- | This contains the actual flowgram for a single read.+data ReadBlock = ReadBlock {+      read_header                :: ReadHeader+    -- The data block+    , flowgram                   :: [Flow]+    , flow_index, bases, quality :: ByteString+    }++instance Show ReadBlock where+    show (ReadBlock h f i b q) =+        show h ++ unlines (map ("     "++) +            ["flowgram:\t"++show f+            , "index:\t"++(concat . intersperse " " . map show . B.unpack) i+            , "bases:\t"++BC.unpack b+            , "quality:\t"++(concat . intersperse " " . map show . B.unpack) q+            , ""+            ])

@@ -51,11 +51,12 @@ unbytes :: Integral a => [Word8] -> a unbytes = Data.List.foldr (\x y -> y*256+x) 0 . map fromIntegral -instance Arbitrary Word8 where-    arbitrary = choose (0,255::Int) >>= return . fromIntegral+-- Conflicts with Bio.Util.TestBase+-- instance Arbitrary Word8 where+--    arbitrary = choose (0,255::Int) >>= return . fromIntegral -prop_bswap :: Word8 -> Word8 -> Word8 -> Word8 -> Bool-prop_bswap a1 a2 a3 a4 = (bswap 4 . decode . BB.pack) [a1,a2,a3,a4] == ((decode . BB.pack) [a4,a3,a2,a1] :: Word32)+-- prop_bswap :: Word8 -> Word8 -> Word8 -> Word8 -> Bool+-- prop_bswap a1 a2 a3 a4 = (bswap 4 . decode . BB.pack) [a1,a2,a3,a4] == ((decode . BB.pack) [a4,a3,a2,a1] :: Word32)  -- basic data types data Header = Header { swap :: Bool,  version, count, reserved :: Word32 }

@@ -14,6 +14,8 @@ lines, mylines :: B.ByteString -> [B.ByteString] lines = case length (B.lines $ B.pack "\n") of 1 -> B.lines                                                0 -> mylines+                                               _ -> error "Unknown implementaion of 'lines'!?"+ 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)

@@ -0,0 +1,132 @@+{-# OPTIONS -fglasgow-exts #-}++module Bio.Util.TestBase where++import Control.Monad (liftM)+import System.CPUTime+import System.Time+import Test.QuickCheck+import System.Random+-- import Data.Char (ord)+import Data.Word+import Data.ByteString.Lazy (pack)++import Bio.Sequence.SeqData++data Test = forall t . Testable t => T String t++newtype Nucleotide = N Char deriving Show+newtype Quality    = Q Word8 deriving Show++fromN :: Nucleotide -> Char+fromN (N c) = c++fromQ :: Quality -> Word8+fromQ (Q c) = c++-- | For testing, variable lengths+newtype EST = E Sequence deriving Show+newtype ESTq = Eq Sequence deriving Show+newtype Protein = P Sequence deriving Show++-- | For benchmarking, fixed lengths+newtype EST_short = ES Sequence deriving Show+newtype EST_long  = EL Sequence deriving Show+newtype EST_set  = ESet [Sequence] deriving Show++-- | Take time (CPU and wall clock) and report it+time :: String -> IO () -> IO ()+time msg action = do+    d1 <- getClockTime+    t1 <- getCPUTime+    action+    t2 <- getCPUTime+    d2 <- getClockTime+    putStrLn $ "\n"++msg++", CPU time: " ++ showT (t2-t1) ++ ", wall clock: "+                 ++ timeDiffToString (diffClockTimes d2 d1)++-- | Print a CPUTime difference+showT :: Integral a => a -> String+showT t = show (fromIntegral t/1e12::Double)++"s"++-- | Shamelessly stolen from FPS+integralRandomR :: (Integral a, RandomGen g) => (a,a) -> g -> (a,g)+integralRandomR  (a,b) g = case randomR (fromIntegral a :: Integer,+                                         fromIntegral b :: Integer) g of+                            (x,g') -> (fromIntegral x, g')++-- | Constrained position generators++genOffset :: Gen Offset+genOffset = do isneg <- arbitrary+               nnoff <- genNonNegOffset+               return $ (if isneg then negate else id) nnoff++genNonNegOffset :: Gen Offset+genNonNegOffset = liftM (subtract 1) genPositiveOffset++genPositiveOffset :: Gen Offset+genPositiveOffset = do scale <- chooseInteger (1, 13)+                       liftM fromIntegral $ chooseInteger (1, 2^scale)+    where chooseInteger :: (Integer, Integer) -> Gen Integer+          chooseInteger = choose+++instance Random Word8 where+    randomR = integralRandomR+    random = randomR (minBound,maxBound)++instance Arbitrary Word8 where+    arbitrary = choose (0,255)+    coarbitrary _ = id++instance Arbitrary Nucleotide where+    arbitrary = elements (map N "aaacccgggtttn")+    coarbitrary _ = id++instance Arbitrary Quality where+    arbitrary = do c <- choose (0,60)+                   return (Q c)+    coarbitrary _ = id++instance Arbitrary ESTq where+    arbitrary = do n <- choose (1,100)+                   s <- vector n+                   q <- vector n+                   return $ Eq $ Seq (fromStr "qctest")+                              (fromStr $ map fromN s) (Just $ pack $ map fromQ q)+    coarbitrary _ = id++instance Arbitrary EST where+    arbitrary = do n <- choose (1,100)+                   s <- vector n+                   return $ E $ Seq (fromStr "qctest")+                              (fromStr $ map fromN s) Nothing+    coarbitrary _ = id++instance Arbitrary Char where+    arbitrary = elements (['A'..'Z']++['a'..'z']++" \t\n\r")+    coarbitrary _ = id++instance Arbitrary EST_short where+    arbitrary = do let n = 200+                   s <- vector n+                   q <- vector n+                   return $ ES $ Seq (fromStr "qctest")+                              (fromStr $ map fromN s) (Just $ pack $ map fromQ q)+    coarbitrary _ = id++instance Arbitrary EST_long where+    arbitrary = do let n = 1000+                   s <- vector n+                   q <- vector n+                   return $ EL $ Seq (fromStr "qctest")+                              (fromStr $ map fromN s) (Just $ pack $ map fromQ q)+    coarbitrary _ = id++-- 1000 ESTs of length 1000+instance Arbitrary EST_set where+    arbitrary = do let n = 1000+                   s <- vector n+                   return (ESet $ map (\(EL x) -> x) s)+    coarbitrary _ = id

@@ -1,5 +1,5 @@ Name:                bio-Version:             0.3.3.4+Version:             0.3.5 License:             LGPL License-file:        LICENSE Author:              Ketil Malde@@ -17,15 +17,16 @@                      phd formats.  Rudimentary support for doing alignments - including                      dynamic adjustment of scores based on sequence quality - and Blast                      output parsing.  Partly implemented single linkage clustering, and-                     multiple alignment.-                     .+		     multiple alignment.  Reading Gene Ontology (GO) annotations (GOA) and+		     definitions\/hierarchy.+		     .                      The Darcs repository is at: <http://malde.org/~ketil/biohaskell/biolib>.-Homepage:            http://malde.org/~ketil/+Homepage:            http://blog.malde.org/index.php/the-haskell-bioinformatics-library/  Tested-With:         GHC==6.8.2 Build-Type:          Simple-Build-Depends:       base>3, QuickCheck<2, binary, tagsoup>=0.4, bytestring,-                     containers, array, parallel, parsec+Build-Depends:       base>3, QuickCheck<2, binary, tagsoup>=0.4, bytestring >= 0.9.1,+                     containers, array, parallel, parsec, random, old-time, mtl -- add fps for ghc 6.4.2; change imports in Bio/Sequence/TwoBit.hs if you want QC 2  -- We omit the debian/ and Test/ files because those are for development, not installation.@@ -33,18 +34,26 @@  Exposed-modules:     Bio.Sequence,                      Bio.Sequence.SeqData,-                     Bio.Sequence.Fasta, Bio.Sequence.TwoBit, Bio.Sequence.Phd,+                     Bio.Sequence.Fasta, Bio.Sequence.FastQ,+		     Bio.Sequence.TwoBit, Bio.Sequence.Phd,                      Bio.Sequence.Entropy, Bio.Sequence.HashWord,                      Bio.Sequence.GOA,                      Bio.Sequence.GeneOntology, 		     Bio.Sequence.KEGG,+		     Bio.Sequence.SFF,                      Bio.Alignment.BlastData, Bio.Alignment.BlastFlat,                      Bio.Alignment.Blast, Bio.Alignment.BlastXML,                      Bio.Alignment.AlignData, Bio.Alignment.Matrices,                      Bio.Alignment.SAlign, Bio.Alignment.AAlign, Bio.Alignment.QAlign-                     Bio.Alignment.Multiple, Bio.Alignment.ACE+                     Bio.Alignment.Multiple, Bio.Alignment.ACE,+                     Bio.Alignment.Soap,                      Bio.Clustering,-                     Bio.Util, Bio.Util.Parsex+                     Bio.Util, Bio.Util.Parsex, Bio.Util.TestBase+		 Bio.Location.Strand, Bio.Location.Position,+		 Bio.Location.ContigLocation, Bio.Location.Location, Bio.Location.LocMap,+		 Bio.Location.OnSeq, Bio.Location.SeqLocation, Bio.Location.SeqLocMap,+		 Bio.GFF3.Escape, Bio.GFF3.Feature, Bio.GFF3.FeatureHier, Bio.GFF3.FeatureHierSequences,+		 Bio.GFF3.SGD  extensions:          CPP, ParallelListComp ghc-options:         -Wall -O2 -fexcess-precision -funbox-strict-fields -auto-all