diff --git a/Bio/Alignment/Blast.hs b/Bio/Alignment/Blast.hs --- a/Bio/Alignment/Blast.hs +++ b/Bio/Alignment/Blast.hs @@ -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" diff --git a/Bio/Alignment/BlastXML.hs b/Bio/Alignment/BlastXML.hs --- a/Bio/Alignment/BlastXML.hs +++ b/Bio/Alignment/BlastXML.hs @@ -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 " 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) diff --git a/Bio/Alignment/Multiple.hs b/Bio/Alignment/Multiple.hs --- a/Bio/Alignment/Multiple.hs +++ b/Bio/Alignment/Multiple.hs @@ -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 diff --git a/Bio/Alignment/Soap.hs b/Bio/Alignment/Soap.hs new file mode 100644 --- /dev/null +++ b/Bio/Alignment/Soap.hs @@ -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) diff --git a/Bio/GFF3/Escape.hs b/Bio/GFF3/Escape.hs new file mode 100644 --- /dev/null +++ b/Bio/GFF3/Escape.hs @@ -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 diff --git a/Bio/GFF3/Feature.hs b/Bio/GFF3/Feature.hs new file mode 100644 --- /dev/null +++ b/Bio/GFF3/Feature.hs @@ -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 diff --git a/Bio/GFF3/FeatureHier.hs b/Bio/GFF3/FeatureHier.hs new file mode 100644 --- /dev/null +++ b/Bio/GFF3/FeatureHier.hs @@ -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 = [] diff --git a/Bio/GFF3/FeatureHierSequences.hs b/Bio/GFF3/FeatureHierSequences.hs new file mode 100644 --- /dev/null +++ b/Bio/GFF3/FeatureHierSequences.hs @@ -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 diff --git a/Bio/GFF3/SGD.hs b/Bio/GFF3/SGD.hs new file mode 100644 --- /dev/null +++ b/Bio/GFF3/SGD.hs @@ -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 diff --git a/Bio/Location/ContigLocation.hs b/Bio/Location/ContigLocation.hs new file mode 100644 --- /dev/null +++ b/Bio/Location/ContigLocation.hs @@ -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) diff --git a/Bio/Location/LocMap.hs b/Bio/Location/LocMap.hs new file mode 100644 --- /dev/null +++ b/Bio/Location/LocMap.hs @@ -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] + diff --git a/Bio/Location/Location.hs b/Bio/Location/Location.hs new file mode 100644 --- /dev/null +++ b/Bio/Location/Location.hs @@ -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 diff --git a/Bio/Location/OnSeq.hs b/Bio/Location/OnSeq.hs new file mode 100644 --- /dev/null +++ b/Bio/Location/OnSeq.hs @@ -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" diff --git a/Bio/Location/Position.hs b/Bio/Location/Position.hs new file mode 100644 --- /dev/null +++ b/Bio/Location/Position.hs @@ -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 -> "(-)" diff --git a/Bio/Location/SeqLocMap.hs b/Bio/Location/SeqLocMap.hs new file mode 100644 --- /dev/null +++ b/Bio/Location/SeqLocMap.hs @@ -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 diff --git a/Bio/Location/SeqLocation.hs b/Bio/Location/SeqLocation.hs new file mode 100644 --- /dev/null +++ b/Bio/Location/SeqLocation.hs @@ -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 + + diff --git a/Bio/Location/Strand.hs b/Bio/Location/Strand.hs new file mode 100644 --- /dev/null +++ b/Bio/Location/Strand.hs @@ -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 + diff --git a/Bio/Sequence.hs b/Bio/Sequence.hs --- a/Bio/Sequence.hs +++ b/Bio/Sequence.hs @@ -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 diff --git a/Bio/Sequence/FastQ.hs b/Bio/Sequence/FastQ.hs new file mode 100644 --- /dev/null +++ b/Bio/Sequence/FastQ.hs @@ -0,0 +1,70 @@ +{- | + Module: Bio.Sequence.FastQ + + Support the FastQ format that combines sequence and quality. See: + + * + + 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: + + * + + * + + 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!") diff --git a/Bio/Sequence/Fasta.hs b/Bio/Sequence/Fasta.hs --- a/Bio/Sequence/Fasta.hs +++ b/Bio/Sequence/Fasta.hs @@ -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. diff --git a/Bio/Sequence/KEGG.hs b/Bio/Sequence/KEGG.hs --- a/Bio/Sequence/KEGG.hs +++ b/Bio/Sequence/KEGG.hs @@ -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) diff --git a/Bio/Sequence/SFF.hs b/Bio/Sequence/SFF.hs new file mode 100644 --- /dev/null +++ b/Bio/Sequence/SFF.hs @@ -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 + , "" + ]) diff --git a/Bio/Sequence/TwoBit.hs b/Bio/Sequence/TwoBit.hs --- a/Bio/Sequence/TwoBit.hs +++ b/Bio/Sequence/TwoBit.hs @@ -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 } diff --git a/Bio/Util.hs b/Bio/Util.hs --- a/Bio/Util.hs +++ b/Bio/Util.hs @@ -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) diff --git a/Bio/Util/TestBase.hs b/Bio/Util/TestBase.hs new file mode 100644 --- /dev/null +++ b/Bio/Util/TestBase.hs @@ -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 diff --git a/bio.cabal b/bio.cabal --- a/bio.cabal +++ b/bio.cabal @@ -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: . -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