bio 0.3.5 → 0.4

raw patch · 27 files changed

+1330/−357 lines, 27 filesdep ~basePVP ok

version bump matches the API change (PVP)

Dependency ranges changed: base

API changes (from Hackage documentation)

Files

• Bio/Alignment/AAlign.hs +6/−6
• Bio/Alignment/ACE.hs +8/−6
• Bio/Alignment/AlignData.hs +5/−5
• Bio/Alignment/Bowtie.hs +153/−0
• Bio/Alignment/Multiple.hs +1/−1
• Bio/Alignment/QAlign.hs +12/−12
• Bio/Alignment/SAlign.hs +8/−8
• Bio/GFF3/FeatureHierSequences.hs +3/−3
• Bio/GFF3/SGD.hs +2/−2
• Bio/Location/ContigLocation.hs +85/−40
• Bio/Location/LocMap.hs +39/−25
• Bio/Location/Location.hs +78/−34
• Bio/Location/OnSeq.hs +38/−9
• Bio/Location/Position.hs +37/−6
• Bio/Location/SeqLocMap.hs +31/−3
• Bio/Location/SeqLocation.hs +43/−5
• Bio/Location/Strand.hs +10/−2
• Bio/Sequence.hs +3/−3
• Bio/Sequence/FastQ.hs +9/−7
• Bio/Sequence/Fasta.hs +19/−15
• Bio/Sequence/Phd.hs +10/−9
• Bio/Sequence/SFF.hs +55/−16
• Bio/Sequence/SFF_name.hs +86/−0
• Bio/Sequence/SeqData.hs +232/−19
• Bio/Sequence/TwoBit.hs +341/−107
• Bio/Util/TestBase.hs +6/−6
• bio.cabal +10/−8

@@ -26,17 +26,17 @@ -- Edit distances  -- | Calculate global edit distance (Needleman-Wunsch alignment score)-global_score :: (Num a, Ord a) => SubstMx a -> (a,a) -> Sequence -> Sequence -> a+global_score :: (Num a, Ord a) => SubstMx t a -> (a,a) -> Sequence t -> Sequence t -> a global_score mx g s1 s2 = uncurry max . last . last                            $ columns (score_select minf mx g) (0,fst g) s1 s2  -- | Calculate local edit distance (Smith-Waterman alignment score)-local_score :: (Num a, Ord a) => SubstMx a -> (a,a) -> Sequence -> Sequence -> a+local_score :: (Num a, Ord a) => SubstMx t a -> (a,a) -> Sequence t -> Sequence t -> a local_score mx g s1 s2 = maximum . map (uncurry max) . concat                           $ columns (score_select 0 mx g) (0,fst g) s1 s2  -- | Generic scoring and selection function for global and local scoring-score_select :: (Num a,Ord a) => a -> SubstMx a -> (a,a) -> Selector (a,a)+score_select :: (Num a,Ord a) => a -> SubstMx t a -> (a,a) -> Selector (a,a) score_select minf mx (go,ge) cds =      let (reps,ids) = partition (isRepl.snd) cds          s = maximum $ minf:[max sub gap +mx (x,y) | ((sub,gap),Repl x y) <- reps]@@ -51,12 +51,12 @@ fp (x,ax) (s,e) = (x+s,e:ax)  -- | Calculate global alignment (Needleman-Wunsch)-global_align :: (Num a, Ord a) => SubstMx a -> (a,a) -> Sequence -> Sequence -> (a,EditList)+global_align :: (Num a, Ord a) => SubstMx t a -> (a,a) -> Sequence t -> Sequence t -> (a,EditList) global_align mx g s1 s2 = revsnd . uncurry max' . last . last                 $ columns (align_select minf mx g) ((0,[]),(fst g,[])) s1 s2  -- | Calculate local alignmnet (Smith-Waterman)-local_align :: (Num a, Ord a) => SubstMx a -> (a,a) -> Sequence -> Sequence -> (a,EditList)+local_align :: (Num a, Ord a) => SubstMx t a -> (a,a) -> Sequence t -> Sequence t -> (a,EditList) local_align mx g s1 s2 = revsnd . maximumBy (compare `on` fst)                          . map (uncurry max') . concat                          $ columns (align_select 0 mx g) ((0,[]),(fst g,[])) s1 s2@@ -64,7 +64,7 @@ revsnd (s,a) = (s,reverse a)  -- | Generic scoring and selection for global and local alignment-align_select :: (Num a, Ord a) => a -> SubstMx a -> (a,a) -> Selector ((a,EditList),(a,EditList))+align_select :: (Num a, Ord a) => a -> SubstMx t a -> (a,a) -> Selector ((a,EditList),(a,EditList)) align_select minf mx (go,ge) cds =      let (reps,ids) = partition (isRepl.snd) cds          s = maximumBy (compare `on` fst) 

@@ -25,6 +25,8 @@           DS (phred header? left empty by CAP3)           RD read2 ...    @++   As far as I know, this is only used for nucleotide sequences. -}  {-# LANGUAGE CPP #-}@@ -47,11 +49,11 @@ import Control.Monad (liftM) -- ,when? import Data.Char (chr) -data Assembly = Asm { contig :: (Sequence,Gaps), fragments :: Alignment }+data Assembly = Asm { contig :: (Sequence Nuc,Gaps), fragments :: Alignment Nuc}                 deriving Show  {-# DEPRECATED reads "Stupid name, replaced by 'fragments'." #-}-reads :: Assembly -> Alignment+reads :: Assembly -> Alignment Nuc reads = fragments -- deprecated, stupid, stupid name.  type Str = ByteString@@ -121,7 +123,7 @@                                 _ -> Nothing) <?> "empty line"  -- | parse the contig and quality information (CO, BQ)-ctg :: AceParser (Sequence,Gaps)+ctg :: AceParser (Sequence Nuc,Gaps) ctg = do   name <- co   sd <- sdata@@ -163,7 +165,7 @@ bq = parse1 (\t -> case t of BQ -> Just (); _ -> Nothing) <?> "BQ"  -- | Given the CO info, get the AFS'es-asm :: (Sequence,Gaps) -> AceParser Assembly+asm :: (Sequence Nuc,Gaps) -> AceParser Assembly asm cg = do   many blank   afs <- many1 af@@ -191,14 +193,14 @@ readInt' :: Str -> Maybe Int readInt' = liftM (fromIntegral . fst) . B.readInt -rds :: (Sequence,Gaps) -> [(Str,Dir,Offset)] -> AceParser Assembly+rds :: (Sequence Nuc,Gaps) -> [(Str,Dir,Offset)] -> AceParser Assembly rds cg xs = do      r <- many1 rseq      -- todo: check the number and merge with the afs      let f (_name,d,off) (s,gs) = (off,d,s,gs)      return $ Asm { contig = cg, fragments = zipWith f xs r } -rseq :: AceParser (Sequence,Gaps)+rseq :: AceParser (Sequence Nuc,Gaps) rseq = do   (rn,_len,_,_) <- rd   (s,gaps) <- return . extractGaps =<< sdata

@@ -37,7 +37,7 @@  data Dir = Fwd | Rev deriving (Eq,Show) type Gaps = [Offset]-type Alignment = [(Offset,Dir,Sequence,Gaps)]+type Alignment a = [(Offset,Dir,Sequence a,Gaps)]  -- | Gaps are coded as '*'s, this function removes them, and returns --   the sequence along with the list of gap positions.@@ -84,11 +84,11 @@ isRepl _ = False  -- | A substitution matrix gives scores for replacing a character with another.---   Typically, it will be symmetric.-type SubstMx a = (Chr,Chr) -> a+--   Typically, it will be symmetric.  It is type-tagged with the alphabet - Nuc or Amino.+type SubstMx t a = (Chr,Chr) -> a  -- | Evaluate an Edit based on SubstMx and gap penalty-eval :: SubstMx a -> a -> Edit -> a+eval :: SubstMx t a -> a -> Edit -> a eval mx g c = case c of Ins _ -> g; Del _ -> g; Repl x y -> mx (x,y)  -- | A Selector consists of a zero element, and a funcition@@ -99,7 +99,7 @@ -- | Calculate a set of columns containing scores --   This represents the columns of the alignment matrix, but will only require linear space --   for score calculation.-columns :: Selector a -> a -> Sequence -> Sequence -> [[a]]+columns :: Selector a -> a -> Sequence b -> Sequence b -> [[a]] columns f z (Seq _ s1 _) (Seq _ s2 _) = columns' f z s1 s2  columns' :: Selector a -> a -> SeqData -> SeqData -> [[a]]

@@ -0,0 +1,153 @@+{-| This module provides a data type to represent an alignment+produced by the Bowtie short-read alignment tool (see+<http://bowtie-bio.sourceforge.net/index.shtml>).++The simple accessors recapitulate the details of the Bowtie alignment+output.  The position of the alignment is given by the \"0-based offset+into the reference sequence where leftmost character of the alignment+occurs\".  Thus, for forward-strand alignments this is the 5\' end of+the query sequence while for reverse-complement alignments this is the+3\' end of the query sequence.  Similarly, the query sequence and query+quality are shown in reference forward strand orientation, and thus+may be reverse complemented.++-}++module Bio.Alignment.Bowtie ( -- * Data type and basic accessors+                              Align(..), Mismatch(..), length, nmismatch, querySequ, queryQual++                            -- * Sequence positions of alignments+                            , refCLoc, refCSeqLoc, refSeqLoc, refSeqPos, mismatchSeqPos++                            -- * Parsing Bowtie output+                            , parse++                            -- * Other utilities+                            , sameRead+                            )+    where ++import Prelude hiding (length)+import Control.Monad.Error+import Control.Monad.State+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.Maybe++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++data Align = Align { name :: !SeqName -- ^ Name of the query sequence+                   , strand :: !Strand -- ^ Strand of the alignment on the reference sequence+                   , refname :: !SeqName -- ^ Name of the reference sequence+                   , leftoffset :: !Offset -- ^ Zero-based offset of the left-most aligned position in the reference+                   , sequ :: !SeqData -- ^ Query sequence, in the reference forward strand orientation+                   , qual :: !QualData -- ^ Query quality, in the reference forward strand orientation+                   , mismatches :: ![Mismatch] -- ^ Mismatches+                   } deriving (Read, Show, Eq, Ord)++-- | Returns the length of the query sequence+length :: Align -> Offset+length = LBS.length . sequ++-- | Returns the number of mismatches in the alignment+nmismatch :: Align -> Int+nmismatch = List.length . mismatches++-- | Parses a line of Bowtie output to produce a 'Align'+parse :: LBS.ByteString -> Either String Align+parse bstr+    = case LBS.split '\t' bstr of+        [nameBStr,strandBStr,refnameBStr,leftoffBStr,sequBStr,qualBStr,_,mismatchesBStr]+            -> do str <- parseStrand strandBStr+                  loff <- parseInt64 leftoffBStr+                  mms <- mapM parseMismatch . LBS.split (',') $ mismatchesBStr+                  return $ Align (LBS.copy nameBStr) str (LBS.copy refnameBStr) loff (LBS.copy sequBStr) (LBS.copy qualBStr) mms+        _ -> throwError $ strMsg $ "Malformed Bowtie alignment " ++ (show . LBS.unpack) bstr+    where parseStrand str | str == (LBS.singleton '+') = return Fwd+                          | str == (LBS.singleton '-') = return RevCompl+                          | otherwise = throwError $ "Unknown strand " ++ (show . LBS.unpack $ str)++-- | Query sequence as given in the query file+querySequ :: Align -> SeqData+querySequ ba = stranded (strand ba) (sequ ba)++-- | Query quality as given in the query file+queryQual :: Align -> QualData+queryQual ba = case strand ba of+                 Fwd -> qual ba+                 RevCompl -> LBSW.reverse $ qual ba++-- | Returns the sequence position of the start of the query sequence+-- alignment.  This will include the strand of the alignment and will+-- not be the same as the position computed from 'leftoffset' when the+-- alignment is on the reverse complement strand.+refSeqPos :: Align -> SeqLoc.SeqPos+refSeqPos ba = OnSeq (refname ba) $ CLoc.startPos $ refCLoc ba++-- | Returns the sequence location covered by the query in+-- the alignment.  This will be a sequence location on the reference+-- sequence and may run on the forward or the reverse complement+-- strand.+refCSeqLoc :: Align -> SeqLoc.ContigSeqLoc+refCSeqLoc ba = OnSeq (refname ba) (refCLoc ba)++-- | As 'refCSeqLoc' but without the reference sequence name.+refCLoc :: Align -> CLoc.ContigLoc+refCLoc ba = CLoc.ContigLoc (leftoffset ba) (length ba) (strand ba)++-- | Returns the sequence location covered by the query, as+-- 'refCSeqLoc', as a 'SeqLoc.SeqLoc' location.+refSeqLoc :: Align -> SeqLoc.SeqLoc +refSeqLoc ba = OnSeq (refname ba) (Loc.Loc [ refCLoc ba ])++-- | Returns true when two alignments were derived from the same+-- sequencing read.  As Bowtie writes alignments of query sequences in+-- their order in the query file, all alignments of a given read are+-- grouped together and the lists of all alignments for each read can+-- be gathered with+-- +-- > groupBy sameRead+sameRead :: Align -> Align -> Bool+sameRead ba1 ba2 = (name ba1) == (name ba2)++-- | Representation of a single mismatch in a bowtie alignment+data Mismatch = Mismatch { mmoffset :: !Offset -- ^ Offset of the mismatch site from the 5\' end of the query+                         , refbase :: !Char -- ^ Reference nucleotide+                         , readbase :: !Char -- ^ Query nucleotide+                         } deriving (Read, Show, Eq, Ord)++-- | Sequence position of a mismatch on the reference sequence.+mismatchSeqPos :: Align -> Mismatch -> SeqLoc.SeqPos+mismatchSeqPos ba mm = OnSeq (refname ba) mmpos+    where mmpos = fromMaybe badOffset $ (Pos.Pos (mmoffset mm) Fwd) `CLoc.posOutof` (refCLoc ba)+          badOffset = error $ "Bad mismatch offset: " ++ show (ba, mm)++parseMismatch :: LBS.ByteString -> Either String Mismatch+parseMismatch = evalStateT parser+    where parser = liftM3 Mismatch parseOffset parseRefNt parseReadNt+          parseOffset = parseTo ':' >>= lift . parseInt64+          parseRefNt = parseTo '>' >>= parseChar+          parseReadNt = parseToEnd >>= parseChar+          parseTo ch = StateT $ \str ->+                       case LBS.span (/= ch) str of+                         (out, rest) -> do (_, afterCh) <- maybe missingChError Right . LBS.uncons $ rest+                                           return (out, afterCh)+              where missingChError = Left ("Failed to find " ++ show ch)+          parseToEnd = StateT $ \str -> return (str, LBS.empty)+          parseChar str = lift $ case LBS.uncons str of+                                   Just (ch, rest) | LBS.null rest -> return ch+                                   _ -> Left ("Malformed character " ++ (show . LBS.unpack) str) ++parseInt64 :: LBS.ByteString -> Either String Offset+parseInt64 zstr = case LBS.readInteger zstr of+                    Just (z, rest) | LBS.null rest -> return $ fromIntegral z+                    _ -> throwError $ "Malformed integer " ++ show zstr+

@@ -14,7 +14,7 @@ -- | Progressive multiple alignment. --   Calculate a tree from agglomerative clustering, then align --   at each branch going bottom up.  Returns a list of columns (rows?).-progressive :: (Sequence -> Sequence -> (Double,EditList)) -> [Sequence] -> [String]+progressive :: (Sequence a -> Sequence a -> (Double,EditList)) -> [Sequence a] -> [String] progressive = undefined  -- |  Derive alignments indirectly, i.e. calculate A|C using alignments A|B and B|C.

@@ -49,10 +49,10 @@ type QSelector a = [(a,Edit,Qual,Qual)] -> a  -- using 21 as default corresponds to blastn default ratio (1,-3)-columns :: QSelector a -> a -> Sequence -> Sequence -> [[a]]+columns :: QSelector a -> a -> Sequence t -> Sequence t -> [[a]] columns sel z s1 s2 = columns' sel (z,r0,c0) s1' s2'     where (s1',s2') = (zup s1, zup s2)-          zup :: Sequence -> [(Chr,Qual)]+          zup :: Sequence t -> [(Chr,Qual)]           zup (Seq _ sd Nothing) = map (\c -> (c,22)) $ B.unpack sd           zup (Seq _ sd (Just qd)) = zip (B.unpack sd) (B.unpack qd) @@ -88,7 +88,7 @@ minf :: Double minf = -100000000 -type QualMx a = Qual -> Qual -> SubstMx a+type QualMx t a = Qual -> Qual -> SubstMx t a  qualMx :: Qual -> Qual -> (Chr,Chr) -> Double qualMx q1 q2 (x,y) = if isN x || isN y then 0.0 else@@ -116,12 +116,12 @@ -- Edit distances  -- | Calculate global edit distance (Needleman-Wunsch alignment score)-global_score :: QualMx Double -> (Double,Double) -> Sequence -> Sequence -> Double+global_score :: QualMx t Double -> (Double,Double) -> Sequence t -> Sequence t -> Double global_score mx g s1 s2 = uncurry max . last . last                           $ columns (score_select minf mx g) (0,fst g) s1 s2  -- | Calculate local edit distance (Smith-Waterman alignment score)-local_score :: QualMx Double -> (Double,Double) -> Sequence -> Sequence -> Double+local_score :: QualMx t Double -> (Double,Double) -> Sequence t -> Sequence t -> Double local_score mx g s1 s2 = maximum . map (uncurry max) . concat                          $ columns (score_select 0 mx g) (0,fst g) s1 s2 @@ -132,7 +132,7 @@ -- Oh. local score_select will not work, since we cannot replace any matrix entry -- with zero in order to initiate a new alignment.  So we need 'minf', except in the -- initial row/column.  Damn.-overlap_score :: QualMx Double -> (Double,Double) -> Sequence -> Sequence -> Double+overlap_score :: QualMx t Double -> (Double,Double) -> Sequence t -> Sequence t -> Double overlap_score mx g s1 s2 = maximum $ map (uncurry max) $ sel cols     where cols   = columns overlap_score_select (0,fst g) s1 s2            sel cs = map last cs ++ last cs@@ -141,7 +141,7 @@           overlap_score_select [_] = (0,minf)  -- | Generic scoring and selection function for global and local scoring-score_select :: Double -> QualMx Double -> (Double,Double) -> QSelector (Double,Double)+score_select :: Double -> QualMx t Double -> (Double,Double) -> QSelector (Double,Double) score_select minf mx (go,ge) cds =     let (reps,ids) = partition (isRepl.snd') cds         s = maximum $ minf:[max sub gap + mx q1 q2 (x,y) | ((sub,gap),Repl x y,q1,q2) <- reps]@@ -161,13 +161,13 @@ fp (x,ax) (s,e) = (x+s,e:ax)  -- | Calculate global alignment (Needleman-Wunsch)-global_align :: QualMx Double -> (Double,Double) -> Sequence -> Sequence -> (Double,EditList)+global_align :: QualMx t Double -> (Double,Double) -> Sequence t -> Sequence t -> (Double,EditList) global_align mx g s1 s2 = revsnd . uncurry max' . last . last                $ columns (align_select minf mx g) ((0,[]),(fst g,[])) s1 s2  -- | Calculate local alignment (Smith-Waterman) --   (can we replace uncurry max' with fst - a local alignment must always end on a subst, no?)-local_align :: QualMx Double -> (Double,Double) -> Sequence -> Sequence -> (Double, EditList)+local_align :: QualMx t Double -> (Double,Double) -> Sequence t -> Sequence t -> (Double, EditList) local_align mx g s1 s2 = revsnd . maximumBy (compare `on` fst)                          . map (uncurry max') . concat                          $ columns (align_select 0 mx g) ((0,[]),(fst g,[])) s1 s2@@ -175,7 +175,7 @@ -- | Calucalte best overlap score, where gaps at the edges are free --   The starting point is like for local score (0 cost for initial indels), --   the result is the maximum anywhere in the last column or bottom row of the matrix.-overlap_align :: QualMx Double -> (Double,Double) -> Sequence -> Sequence -> (Double,EditList)+overlap_align :: QualMx t Double -> (Double,Double) -> Sequence t -> Sequence t -> (Double,EditList) overlap_align mx g s1 s2 = revsnd . maximumBy (compare `on` fst) . map (uncurry max') $ sel cols     where cols   = columns overlap_align_select ((0,[]),(minf,[])) s1 s2            sel cs = map last cs ++ last cs@@ -184,7 +184,7 @@           overlap_align_select [_] = ((0,[]),(minf,[]))  -- | Variant that retains indels to retain the entire sequence in the result-overlap_align' :: QualMx Double -> (Double,Double) -> Sequence -> Sequence -> (Double,EditList)+overlap_align' :: QualMx t Double -> (Double,Double) -> Sequence t -> Sequence t -> (Double,EditList) overlap_align' mx g s1 s2 = revsnd . maximumBy (compare `on` fst) . map (uncurry max') $ sel cols     where cols   = columns overlap_align_select ((0,[]),(fst g,[])) s1 s2            sel cs = map last cs ++ last cs@@ -196,7 +196,7 @@ revsnd (s,a) = (s,reverse a)  -- | Generic scoring and selection for global and local alignment-align_select :: Double -> QualMx Double -> (Double,Double) -> QSelector ((Double,EditList),(Double,EditList))+align_select :: Double -> QualMx t Double -> (Double,Double) -> QSelector ((Double,EditList),(Double,EditList)) align_select minf mx (go,ge) cds =     let (reps,ids) = partition (isRepl.snd') cds         s = maximumBy (compare `on` fst)

@@ -22,19 +22,19 @@ -- Edit distances  -- | Calculate global edit distance (Needleman-Wunsch alignment score)-global_score :: (Num a, Ord a) => SubstMx a -> a -> Sequence -> Sequence -> a+global_score :: (Num a, Ord a) => SubstMx t a -> a -> Sequence t -> Sequence t -> a global_score mx g s1 s2 = last . last $ columns (g_score mx g) 0 s1 s2  -- | Scoring\/selection function for global alignment-g_score:: (Num a,Ord a) => SubstMx a -> a -> Selector a+g_score:: (Num a,Ord a) => SubstMx t a -> a -> Selector a g_score mx g cds = maximum [s+eval mx g e | (s,e) <- cds]  -- | Calculate local edit distance (Smith-Waterman alignment score)-local_score :: (Num a, Ord a) => SubstMx a -> a -> Sequence -> Sequence -> a+local_score :: (Num a, Ord a) => SubstMx t a -> a -> Sequence t -> Sequence t -> a local_score mx g s1 s2 = maximum . concat $ columns (l_score mx g) 0 s1 s2  -- | Scoring\/selection funciton for local alignmnet-l_score :: (Num a,Ord a) => SubstMx a -> a -> Selector a+l_score :: (Num a,Ord a) => SubstMx t a -> a -> Selector a l_score mx g cds = maximum (0:[s+eval mx g e | (s,e) <- cds])  -- ------------------------------------------------------------@@ -46,17 +46,17 @@ -- the score in each cell.  Unreachable cells can then be GC'ed.  -- | Calculate alignments.-global_align :: (Num a, Ord a) => SubstMx a -> a -> Sequence -> Sequence -> EditList+global_align :: (Num a, Ord a) => SubstMx t a -> a -> Sequence t -> Sequence t -> EditList global_align mx g s1 s2 = reverse . snd . last . last                            $ columns (g_align mx g) (0,[]) s1 s2 -g_align :: (Num a, Ord a) => SubstMx a -> a -> Selector (a,EditList)+g_align :: (Num a, Ord a) => SubstMx t a -> a -> Selector (a,EditList) g_align mx g cds = maximumBy (compare `on` fst) [(s+eval mx g e,e:a) | ((s,a),e) <- cds] -local_align :: (Num a, Ord a) => SubstMx a -> a -> Sequence -> Sequence -> EditList +local_align :: (Num a, Ord a) => SubstMx t a -> a -> Sequence t -> Sequence t -> EditList  local_align mx g s1 s2 = reverse . snd . maximumBy (compare `on` fst) . concat                           $ columns (l_align mx g) (0,[]) s1 s2 -l_align :: (Num a, Ord a) => SubstMx a -> a -> Selector (a,EditList)+l_align :: (Num a, Ord a) => SubstMx t a -> a -> Selector (a,EditList) l_align mx g cds = maximumBy (compare `on` fst) $ (0,[]):[(s+eval mx g e,e:a) | ((s,a),e) <- cds] 

@@ -32,7 +32,7 @@ parse str = do (feats, seqlines) <- F.parseWithFasta str                fromLists feats $ mkSeqs seqlines -fromLists :: (Error e, MonadError e m) => [F.Feature] -> [Sequence] -> m FeatureHierSequences+fromLists :: (Error e, MonadError e m) => [F.Feature] -> [Sequence a] -> 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@@ -44,7 +44,7 @@ features :: FeatureHierSequences -> S.Set F.Feature features = FH.features . featureHier -sequences :: FeatureHierSequences -> [Sequence]+sequences :: FeatureHierSequences -> [Sequence a] sequences = M.foldWithKey mkSeqList [] . sequenceMap     where mkSeqList seqname sequ = ((Seq seqname sequ Nothing) :) @@ -65,7 +65,7 @@ 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 :: (Error e, MonadError e m) => FeatureHierSequences -> F.Feature -> m (Sequence a) 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

@@ -46,7 +46,7 @@     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 :: (Error e, MonadError e m) => FeatureHierSequences -> F.Feature -> m (Sequence a) 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@@ -58,7 +58,7 @@ 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 :: (Error e, MonadError e m) => FeatureHierSequences -> F.Feature -> m (Sequence a) 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

@@ -1,15 +1,29 @@--- |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'.+{-| Data type for a sequence location consiting of a contiguous range+of positions on the sequence. -module Bio.Location.ContigLocation ( ContigLoc(..), fromStartEnd, fromPosLen-                                   , bounds, startPos, endPos-                                   , slide, extend, posInto, posOutof-                                   , seqData, seqDataPadded, isWithin, overlaps-                                   , display-                                   )-    where +Throughout, /sequence position/ refers to a 'Pos.Pos' which includes a+strand.  An index into a sequence is referred to as an /offset/, and+is generally of type 'Offset'. + -}++module Bio.Location.ContigLocation ( +  -- * Sequence locations+  ContigLoc(..), fromStartEnd, fromPosLen++  -- * Locations and positions+  , bounds, startPos, endPos, posInto, posOutof, isWithin, overlaps++  -- * Extracting subsequences+  , seqData, seqDataPadded++  -- * Transforming locations+  , slide, extend++  -- * Displaying locations+  , display+                                   ) where+ import Prelude hiding (length)  import Control.Monad.Error@@ -20,56 +34,69 @@ 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+-- | Contiguous sequence location defined by a span of sequence+-- positions, lying on a specific strand of the sequence.+data ContigLoc = ContigLoc { offset5 :: !Offset   -- ^ The offset of the 5\' end of the location, as a 0-based index+                           , length :: !Offset    -- ^ The length of the location+                           , strand :: !Strand    -- ^ The strand of the location                            } 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.+-- | Create a sequence location lying between 0-based starting and+-- ending offsets.  When @start < end@, the location+-- be on the forward strand, otherwise it will be on the+-- reverse complement 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.+-- | Create a sequence location from the sequence position of the+-- start of the location and the length of the position.  The strand+-- of the location, and the direction it extends from the starting+-- position, are determined by the strand of the starting position. 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.+-- | The bounds of a sequence location.  This is a pair consisting of+-- the lowest and highest sequence offsets covered by the region.  The+-- bounds ignore the strand of the sequence location, and 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+-- | Sequence position of the start of the location.  This is the 5'+-- end on the location strand, which will have a higher offset than+-- 'endPos' if the location is on the 'RevCompl' strand. 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+-- | Sequence position of the end of the location, as described in 'startPos'. 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+-- | Returns a location resulting from sliding the original location+-- along the sequence by a specified offset.  A positive offset will+-- move the location away from the 5\' end of the forward stand of the+-- sequence regardless of the strand of the location itself.  Thus,+-- +-- > slide (revCompl cloc) off == revCompl (slide cloc off) 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.+-- | Extract the nucleotide 'SeqData' for the sequence location.  If+-- any part of the location lies outside the bounds of the sequence,+-- an error results. seqData :: (Error e, MonadError e m) => SeqData -> ContigLoc -> m SeqData seqData sequ (ContigLoc seq5 len str)     | seq5 < 0 = outOfBounds@@ -78,7 +105,9 @@                            | 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+-- | As 'seqData', extract the nucleotide subsequence for the+-- location.  Any positions in the location lying outside the bounds+-- of the sequence are returned as @N@ rather than producing an error. seqDataPadded :: SeqData -> ContigLoc -> SeqData seqDataPadded sequ (ContigLoc seq5 len str) = stranded str fwdseq     where fwdseq@@ -90,9 +119,12 @@               | 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'.+-- | Given a sequence position and a sequence location relative to the+-- same sequence, compute a new position representing the original+-- position relative to the subsequence defined by the location.  If+-- the sequence position lies outside of the sequence location,+-- @Nothing@ is returned; thus, the offset of the new position will+-- always be in the range @[0, length cloc - 1]@. posInto :: Pos.Pos -> ContigLoc -> Maybe Pos.Pos posInto (Pos.Pos pos pStrand) (ContigLoc seq5 len cStrand)     | pos < seq5 || pos >= seq5 + len = Nothing@@ -100,9 +132,14 @@                            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'.++-- | Given a sequence location and a sequence position within that+-- location, compute a new position representing the original position+-- relative to the outer sequence.  If the sequence position lies+-- outside the location, @Nothing@ is returned.+-- +-- This function inverts 'posInto' when the sequence position lies+-- within the position is actually within the location. posOutof :: Pos.Pos -> ContigLoc -> Maybe Pos.Pos posOutof (Pos.Pos pos pStrand) (ContigLoc seq5 len cStrand)     | pos < 0 || pos >= len = Nothing@@ -110,26 +147,34 @@                            Fwd      -> Pos.Pos (pos + seq5)             pStrand                            RevCompl -> Pos.Pos (seq5 + len - (pos + 1)) (revCompl pStrand) --- | 'ContigLoc' extended on the 5' and 3' ends.+-- | Returns a sequence location produced by extending the original+-- location on each end, based on a pair of (/5\' extension/, /3\'+-- extension/).  The 5\' extension is applied to the 5\' end of the+-- location on the location strand; if the location is on the+-- 'RevCompl' strand, the 5\' end will have a higher offset than the+-- 3\' end and this offset will increase by the amount of the 5\'+-- extension.  Similarly, the 3\' extension is applied to the 3\' end+-- of the location. 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'.+-- | Returns @True@ when  a sequence position lies within a sequence+-- location on the same sequence, and occupies the same strand. 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.+-- | Returns @True@ when two sequence locations overlap at any+-- position. 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 a human-friendly, zero-based representation of a sequence location. display :: ContigLoc -> String display cloc = show (Pos.offset $ startPos cloc) ++ "to" ++ show (Pos.offset $ endPos cloc)

@@ -1,12 +1,28 @@--- | 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 +{-| Efficient lookup of sequence positions and locations in a large+map of target locations.  For example, target locations might+represent a collection of genes annotated on a chromosome.  The+'LocMap' would efficiently find which gene(s) overlapped a sequence+position on that chromosome. +Target locations are assigned to one or more zones based on+'Loc.bounds'.  Query locations are then tested only against the target+locations in the relevant zones.+-}++module Bio.Location.LocMap (+  -- * Constructing location lookup maps+  LocMap, fromList++  -- * Searching for target locations+  , lookupWithin, lookupOverlaps++  -- * Modifying location lookup maps+  , delete, deleteBy, insert++  -- * Utilities+  , checkInvariants+  ) where + import qualified Data.IntMap as IM import qualified Data.IntSet as IS import Data.List hiding (insert, delete, deleteBy)@@ -20,23 +36,20 @@ defaultZonesize :: Offset defaultZonesize = 65536 --- | Collection mapping a collection of 'Loc' locations, possibly---   overlapping, binned for efficient lookup by position.+-- | Data structure allowing efficient lookup of target sequence+-- locations that overlap a query location.  Target locations can be+-- paired with an arbitrary object. data LocMap a = LocMap !Offset !(IM.IntMap (Loc.Loc, a)) !(IM.IntMap IS.IntSet)  instance Monoid (LocMap a) where-    mempty = mkLocMap defaultZonesize+    mempty = LocMap defaultZonesize IM.empty IM.empty     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.+-- | Create a 'LocMap' from an association list of target locations. fromList :: Offset -> [(Loc.Loc, a)] -> LocMap a-fromList zonesize = foldl' (\lm0 (l,x) -> insert l x lm0) (mkLocMap zonesize)+fromList zonesize = foldl' (\lm0 (l,x) -> insert l x lm0) mempty --- | Add an object with an associated 'Loc' sequence region+-- | Insert a new target association into a target location map. insert :: Loc.Loc -> a -> LocMap a -> LocMap a insert l x (LocMap zonesize keyToLoc zoneKeys)      = let !key = case IM.maxViewWithKey keyToLoc of@@ -50,28 +63,29 @@                                                 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'+-- | Find the (possibly empty) list of target locations and associated+-- objects that contain a sequence position, in the sense of+-- 'Loc.isWithin' 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'+-- | Find the (possibly empty) list of target locations and associated+-- objects that overlap a sequence location, in the sense of 'Loc.overlaps' 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+-- | Remove a target location and 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.+-- | Generalized version of 'delete' that removes the first target+-- location / object association that satisfies 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

@@ -1,21 +1,35 @@--- |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'.+{-| Data type for a more general sequence location consiting of+potentially disjoint ranges of positions on the sequence. -module Bio.Location.Location ( Loc(..), bounds, length, startPos, endPos-                             , extend, posInto, posOutof, isWithin, overlaps-                             , seqData, seqDataPadded-                             , display-                             )-    where +Throughout, /sequence position/ refers to a 'Pos.Pos' which includes a+strand.  An index into a sequence is referred to as an /offset/, and+is generally of type 'Offset'. + -}++module Bio.Location.Location ( +  -- * Sequence locations+  Loc(..)++  -- * Locations and positions+  , bounds, length, startPos, endPos, posInto, posOutof, isWithin, overlaps++  -- * Extracting subsequences+  , seqData, seqDataPadded++  -- * Transforming locations+  , extend++  -- * Displaying locations+  , 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@@ -23,49 +37,65 @@ import Bio.Sequence.SeqData  -- | General (disjoint) sequence region consisting of a concatenated---   set of contiguous regions+--   set of contiguous regions (see 'CLoc.ContigLoc'). 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+-- | Returns the 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.+-- | The bounds of a sequence location.  This is a pair consisting of+-- the lowest and highest sequence offsets covered by the region.  The+-- bounds ignore the strand of the sequence location, and the first+-- element of the pair will always be lower than the second.  Even if+-- the positions in the location do not run monotonically through the+-- location, the overall lowest and highest sequence offsets are returned. 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.+-- | Sequence position of the start of the location.  This is the 5'+-- end on the location strand, which will have a higher offset than+-- 'endPos' if the location is on the 'RevCompl' strand. 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.+-- | Sequence position of the end of the location, as described in 'startPos'. 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.+-- | Extract the nucleotide 'SeqData' for the sequence location.  If+-- any part of the location lies outside the bounds of the sequence,+-- an error results. seqData :: (Error e, MonadError e m) => SeqData -> Loc -> m SeqData seqData sequ (Loc contigs)     = liftM LBS.concat $ mapM (CLoc.seqData sequ) contigs +-- | As 'seqData', extract the nucleotide subsequence for the+-- location.  Any positions in the location lying outside the bounds+-- of the sequence are returned as @N@ rather than producing an error. 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.+-- | Given a sequence position and a sequence location relative to the+-- same sequence, compute a new position representing the original+-- position relative to the subsequence defined by the location.  If+-- the sequence position lies outside of the sequence location,+-- @Nothing@ is returned; thus, the offset of the new position will+-- always be in the range @[0, length cloc - 1]@.+--+-- When the sequence positions in the location are not monotonic,+-- there may be multiple possible posInto solutions.  That is, if the+-- same outer sequence position is covered by two different contiguous+-- blocks of the location, then it would have two possible sequence+-- positions relative to the location. In this case, the position+-- 5\'-most in the location orientation is returned. posInto :: Pos.Pos -> Loc -> Maybe Pos.Pos posInto seqpos (Loc contigs) = posIntoContigs seqpos contigs @@ -76,9 +106,15 @@         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.+-- | Given a sequence location and a sequence position within that+-- location, compute a new position representing the original position+-- relative to the outer sequence.  If the sequence position lies+-- outside the location, @Nothing@ is returned.  +-- +-- This function inverts 'posInto' when the sequence position lies+-- within the position is actually within the location.  Due to the+-- possibility of redundant location-relative positions for a given+-- absolute position, 'posInto' does not necessary invert 'posOutof' posOutof :: Pos.Pos -> Loc -> Maybe Pos.Pos posOutof pos (Loc contigs) = posOutofContigs pos contigs @@ -89,8 +125,15 @@         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+-- | Returns a sequence location produced by extending the original+-- location on each end, based on a pair of (/5\' extension/, /3\'+-- extension/).  These add contiguous positions to the 5\' and 3\'+-- ends of the original location.  The 5\' extension is applied to the+-- 5\' end of the location on the location strand; if the location is+-- on the 'RevCompl' strand, the 5\' end will have a higher offset+-- than the 3\' end and this offset will increase by the amount of the+-- 5\' extension.  Similarly, the 3\' extension is applied to the 3\'+-- end of the location. extend :: (Offset, Offset) -- ^ (5' extension, 3' extension)        -> Loc -> Loc extend _ (Loc []) = error "extendLoc on zero-contig Loc"@@ -101,8 +144,8 @@           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?+-- | Returns @True@ when  a sequence position lies within a sequence+-- location on the same sequence, and occupies the same strand. isWithin :: Pos.Pos -> Loc -> Bool isWithin seqpos (Loc contigs) = or $ map (CLoc.isWithin seqpos) contigs @@ -110,10 +153,11 @@ 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?+-- | Returns @True@ when two sequence locations overlap at any+-- position. overlaps :: Loc -> Loc -> Bool overlaps l1 l2 = not $ null $ overlappingContigs l1 l2 +-- | Display a human-friendly, zero-based representation of a sequence location. display :: Loc -> String display (Loc contigs) = intercalate ";" $ map CLoc.display contigs

@@ -1,9 +1,21 @@-module Bio.Location.OnSeq ( SeqName -                          , OnSeq(..), withSeqData, andSameSeq, onSameSeq-                          , OnSeqs, perSeq, perSeqUpdate, withNameAndSeq-                          )-    where +{-| Data types for functorially lifting sequence positions and+locations onto named sequences.  These are useful for taking functions+that work with sequence positions and locations and associating them+specific, named sequences. +-}++module Bio.Location.OnSeq ( +  -- * Data types+  SeqName , OnSeq(..)++  -- * Utility functions+  , withSeqData, andSameSeq, onSameSeq++  -- * Sequence collections indexed by name+  , OnSeqs, perSeq, perSeqUpdate, withNameAndSeq+  ) where + import Control.Monad.Error import qualified Data.ByteString.Lazy.Char8 as LBS import Data.List@@ -12,8 +24,10 @@  import Bio.Sequence.SeqData +-- | Sequence name, as in a 'Sequence' type SeqName = SeqData +-- | Data type for an object associated with a specific, named sequence data OnSeq a = OnSeq { onSeqName :: !SeqName                      , onSeqObj :: !a                      } deriving (Eq, Ord, Show)@@ -21,27 +35,42 @@ 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+-- | Looks up a sequence by name and applies a function to it+withSeqData :: (Monad m) =>+               (SeqData -> a -> m b) -- ^ Function using sequence data+               -> (SeqName -> m SeqData) -- ^ Lookup sequence by name +               -> OnSeq a -- ^ Object with named sequence +               -> m b withSeqData f lookupSeq (OnSeq seqname x) = lookupSeq seqname >>= flip f x +-- | Tests a predicate when two objects are on the same sequence,+-- returning @False@ if they are on different sequences. 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+-- | Performs an action when two objects are on the same sequence and+-- produces an error otherwise.+onSameSeq :: (Error e, MonadError e 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)-+    | otherwise      = throwError $ strMsg $ "onSameSeq: " ++ show (LBS.unpack xname) ++ " /= " ++ show (LBS.unpack yname) +-- | Data type for a collection of objects indexed by sequence name type OnSeqs a = M.Map SeqName a +-- | Lifts a function on an underlying object to look up the sequence+-- name in a name-indexed collection. perSeq :: (Monoid b) => (a -> b -> c) -> OnSeq a -> OnSeqs b -> c perSeq f (OnSeq seqname x) = f x . M.findWithDefault mempty seqname +-- | Lifts a function that updates an underlying object to look up the+-- named sequence and update a named-index collection. 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 +-- | Lifts a function on underlying objects to look up a sequence in a+-- name-indexed collection 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

@@ -1,8 +1,21 @@--- | Positions on a sequence.  Zero-based 'Int64' indices are used--- throughout, to facilitate direct use of indexing functions on--- 'SeqData'.+{-| Data type for a sequence position. -module Bio.Location.Position ( Pos(..), slide, seqNt, seqNtPadded, display )+Zero-based 'Offset' / 'Int64' indices are used throughout, to+facilitate direct use of indexing functions on 'SeqData'.+-}++module Bio.Location.Position ( +  -- * Sequence positions+  Pos(..)++  -- * Manipulating positions+  , slide++  -- * Extracting sequences+  , seqNt, seqNtPadded++  -- * Displaying positions+  , display )     where   import Control.Monad.Error@@ -15,25 +28,43 @@  -- | Position in a sequence data Pos = Pos { offset :: !Offset -- ^ 0-based index of the position-               , strand :: !Strand -- ^ Optional strand of the position+               , strand :: !Strand -- ^ 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+-- | Returns a position resulting from sliding the original position+-- along the sequence by a specified offset.  A positive offset will+-- move the position away from the 5\' end of the forward stand of the+-- sequence regardless of the strand of the position itself.  Thus,+-- +-- > slide (revCompl pos) off == revCompl (slide pos off)+ slide :: Pos -> Offset -> Pos slide (Pos off str) doff = Pos (off + doff) str +++-- | Extract the nucleotide at a specific sequence position.  If the+-- position lies outside the bounds of the sequence, an error results. 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" +-- | As 'seqNt', extract the nucleotide at a specific sequence+-- position, but return @N@ when the position lies outside the+-- bounds of the sequence.+-- +-- > seqNtPadded sequ pos == (either 'N' id . seqNt sequ) pos seqNtPadded :: SeqData -> Pos -> Char seqNtPadded sequ (Pos off str) | off >= 0 && off < LBS.length sequ = stranded str $ sequ `LBS.index` off                                | otherwise = 'N' +++-- | Display a human-friendly, zero-based representation of a sequence position. display :: Pos -> String display (Pos off str) = show off ++ displayStrand     where displayStrand = case str of

@@ -1,6 +1,19 @@-module Bio.Location.SeqLocMap ( SeqLocMap-                              , empty-                              , fromList, insert+{-| ++Efficient lookup of query positions in a collection of target sequence+locations where positions and locations are associated with specific+sequence names.  This is an extension of 'LocMap' to use locations and+positions on named sequences as in 'SeqLocation'.++-}++module Bio.Location.SeqLocMap ( -- * Location lookup maps for named sequences+                                SeqLocMap, empty, fromList++                              -- * Modifying location lookup maps+                              , insert+                              +                              -- * Searching for target locations                               , lookupWithin, lookupOverlaps                               )     where @@ -13,21 +26,36 @@ import qualified Bio.Location.LocMap as LM import qualified Bio.Location.SeqLocation as SeqLoc +-- | A data structure for efficiently finding target sequence+-- locations ('SeqLoc.Loc') that overlap query positions or locations.+-- Each target location can be associated with an arbitrary additional+-- value in the lookup map. type SeqLocMap a = OnSeqs (LM.LocMap a) +-- | Empty lookup map. empty :: SeqLocMap a empty = M.empty +-- | Creates a 'SeqLocMap' from a list of target locations and their+-- associated objects fromList :: [(SeqLoc.SeqLoc, a)] -> SeqLocMap a fromList = foldl' (\lm (sl,x) -> insert sl x lm) M.empty +-- | Inserts a new target location and associated object into the+-- location lookup map. insert :: SeqLoc.SeqLoc -> a -> SeqLocMap a -> SeqLocMap a insert sloc x = perSeqUpdate (\loc locmap -> LM.insert loc x locmap) sloc +-- | Find the (possibly empty) list of target locations and associated+-- objects that contain a sequence position, in the sense of+-- 'Loc.isWithin'. lookupWithin :: SeqLoc.SeqPos -> SeqLocMap a -> [(SeqLoc.SeqLoc, a)] lookupWithin = withNameAndSeq namedLookupWithin     where namedLookupWithin seqname pos = map (first $ OnSeq seqname) . LM.lookupWithin pos +-- | Find the (possibly empty) list of target locations and associated+-- objects that overlap a sequence location, in the sense of+-- 'Loc.overlaps'. lookupOverlaps :: SeqLoc.SeqLoc -> SeqLocMap a -> [(SeqLoc.SeqLoc, a)] lookupOverlaps = withNameAndSeq namedLookupOverlaps     where namedLookupOverlaps seqname loc = map (first $ OnSeq seqname) . LM.lookupOverlaps loc

@@ -1,14 +1,32 @@-module Bio.Location.SeqLocation ( SeqPos, displaySeqPos-                                , ContigSeqLoc, withinContigSeqLoc, displayContigSeqLoc +{-| ++Data types for sequence locations and sequence positions associated+with specific, named sequences.++-}++module Bio.Location.SeqLocation ( -- * Positions on named sequences+                                  SeqPos++                                -- * Contiguous location spans on named sequences+                                , ContigSeqLoc, withinContigSeqLoc++                                -- * Arbitrary location spans on named sequences                                 , SeqLoc-                                , isWithin, overlaps, seqData-                                , display++                                -- * Testing for location intersection on named sequences+                                , isWithin, overlaps++                                -- * Extracting subsequences+                                , seqData++                                -- * Displaying locations on named sequences+                                , displaySeqPos, displayContigSeqLoc, 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@@ -16,30 +34,50 @@ import qualified Bio.Location.Position as Pos import Bio.Sequence.SeqData +-- | A position on a named sequence type SeqPos = OnSeq Pos.Pos +-- | Display a human-friendly representation of a 'SeqPos' displaySeqPos :: SeqPos -> String displaySeqPos (OnSeq refname pos) = LBS.unpack refname ++ "@" ++ Pos.display pos +-- | A location consisting of a contiguous span of positions on a+-- named sequence. type ContigSeqLoc = OnSeq CLoc.ContigLoc +-- | Test whether a sequence position lies within a sequence location.+-- This requires that the position lie within the location as per+-- 'CLoc.isWithin' and have the same sequence name. withinContigSeqLoc :: SeqPos -> ContigSeqLoc -> Bool withinContigSeqLoc = andSameSeq CLoc.isWithin +-- | Display a human-friendly representation of a 'ContigSeqLoc' displayContigSeqLoc :: ContigSeqLoc -> String displayContigSeqLoc (OnSeq refname cloc) = LBS.unpack refname ++ "@" ++ CLoc.display cloc +-- | A general location, consisting of spans of sequence positions on+-- a specific, named sequence. type SeqLoc = OnSeq Loc.Loc +-- | Test whether a sequence position lies within a sequence location.+-- This requires that the position lie within the location as per+-- 'Loc.isWithin' and have the same sequence name. isWithin :: SeqPos -> SeqLoc -> Bool isWithin = andSameSeq Loc.isWithin +-- | Test whether two sequence locations overlap in any position.+-- This requires that the locations overlap as per 'Loc.overlaps' and+-- have the same sequence name. overlaps :: SeqLoc -> SeqLoc -> Bool overlaps = andSameSeq Loc.overlaps +-- | Extract the subsequence specified by a sequence location from a+-- sequence database.  The sequence name is used to retrieve the full+-- sequence and the subsequence is extracted as by 'Loc.seqData' seqData :: (Error e, MonadError e m) => (SeqName -> m SeqData) -> SeqLoc -> m SeqData seqData = withSeqData Loc.seqData +-- | Display a human-friendly representation of a 'SeqLoc' display :: SeqLoc -> String display (OnSeq refname loc) = LBS.unpack refname ++ "@" ++ Loc.display loc 

@@ -1,3 +1,9 @@+{-| Utilities for manipulating nucleotide sequences and locations on+nucleotide sequences that occur on a forward or a reverse-complement+strand.++-}+ module Bio.Location.Strand ( Strand(..)                            , Stranded(..), stranded                            )@@ -11,11 +17,13 @@ -- | 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.+-- | A nucleotide sequence or location on a nucleotide sequence that+--   lies on a specific strand and has an orientation. class Stranded s where     revCompl :: s -> s +-- | Convert the orientation of a 'Stranded' thing based on a+--   specified 'Strand' stranded :: (Stranded s) => Strand -> s -> s stranded Fwd      = id stranded RevCompl = revCompl

@@ -7,7 +7,7 @@ module Bio.Sequence      (     -- * Data structures etc ("Bio.Sequence.SeqData")-      Sequence(..), Offset, SeqData, Qual, QualData+      Sequence(..), Unknown, Offset, SeqData, Qual, QualData     -- ** Accessor functions     , seqlength, seqlabel, seqheader, seqdata, seqqual, (!)     , appendHeader, setHeader@@ -15,9 +15,9 @@     -- ** Converting to and from String.     , fromStr, toStr     -- ** Nucleotide functionality.-    , compl, revcompl+    , compl, revcompl, revcompl', Nuc, castToNuc     -- ** Protein sequence functionality-    , Amino(..), translate, fromIUPAC, toIUPAC+    , Amino(..), translate, fromIUPAC, toIUPAC, castToAmino      -- * File formats     -- ** The Fasta file format ("Bio.Sequence.Fasta")

@@ -15,6 +15,8 @@     Currently, we only support the non-Solexa FastQ, adding\/subtracting 33 for     the quality values.++   As far as I know, FastQ is only used for nucleotide sequences, never amino acid. -}  module Bio.Sequence.FastQ @@ -32,18 +34,18 @@  import Bio.Sequence.SeqData -readFastQ :: FilePath -> IO [Sequence]+readFastQ :: FilePath -> IO [Sequence Nuc] readFastQ f = (go . B.lines) `fmap` B.readFile f  -hReadFastQ :: Handle -> IO [Sequence]+hReadFastQ :: Handle -> IO [Sequence Nuc] hReadFastQ h = (go . B.lines) `fmap` B.hGetContents h -go :: [B.ByteString] -> [Sequence]+go :: [B.ByteString] -> [Sequence Nuc] 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 :: [B.ByteString] -> Maybe (Either String (Sequence Nuc), [B.ByteString]) parse (h1:sd:h2:sq:rest) =      case (B.uncons h1,B.uncons h2) of       (Just ('@',h1name), Just ('+',h2name))@@ -56,14 +58,14 @@                err = Left $ "Bio.Sequence.FastQ: illegal number of lines in FastQ format: " ++ showStanza            in Just (err, []) -writeFastQ :: FilePath -> [Sequence] -> IO ()+writeFastQ :: FilePath -> [Sequence a] -> IO () writeFastQ f = B.writeFile f . B.concat . map unparse -hWriteFastQ :: Handle -> [Sequence] -> IO ()+hWriteFastQ :: Handle -> [Sequence a] -> IO () hWriteFastQ h = B.hPut h . B.concat . map unparse  -- helper function for writing-unparse :: Sequence -> B.ByteString+unparse :: Sequence a -> 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) 

@@ -5,6 +5,10 @@    sequence data in the Fasta format.    Each sequence consists of a header (with a '>' prefix)    and a set of lines containing the sequence data.++   As Fasta is used for both amino acids and nucleotides, the+   resulting 'Sequence's are type-tagged with 'Unknown'.  If you know the +   type of sequence you are reading, use 'castToAmino' or 'castToNuc'. -}  module Bio.Sequence.Fasta@@ -57,23 +61,23 @@ -}  -- | Lazily read sequences from a FASTA-formatted file-readFasta :: FilePath -> IO [Sequence]+readFasta :: FilePath -> IO [Sequence Unknown] readFasta f = (mkSeqs . B.lines) `fmap` B.readFile f  -- | Write sequences to a FASTA-formatted file. --   Line length is 60.-writeFasta :: FilePath -> [Sequence] -> IO ()+writeFasta :: FilePath -> [Sequence a] -> IO () writeFasta f ss = do   h <- openFile f WriteMode   hWriteFasta h ss   hClose h  -- | Read quality data for sequences to a file.-readQual :: FilePath -> IO [Sequence]+readQual :: FilePath -> IO [Sequence Unknown] readQual f = (mkQual . B.lines) `fmap` B.readFile f  -- | Write quality data for sequences to a file.-writeQual :: FilePath -> [Sequence] -> IO ()+writeQual :: FilePath -> [Sequence a] -> IO () writeQual f ss = do   h <- openFile f WriteMode   hWriteQual h ss@@ -81,7 +85,7 @@  -- | Read sequence and associated quality.  Will error if --   the sequences and qualites do not match one-to-one in sequence.-readFastaQual :: FilePath -> FilePath -> IO [Sequence]+readFastaQual :: FilePath -> FilePath -> IO [Sequence Unknown] readFastaQual  s q = do   ss <- readFasta s   qs <- readQual q@@ -100,7 +104,7 @@  -- | Write sequence and quality data simulatnously --   This may be more laziness-friendly.-writeFastaQual :: FilePath -> FilePath -> [Sequence] -> IO ()+writeFastaQual :: FilePath -> FilePath -> [Sequence a] -> IO () writeFastaQual f q ss = do   hf <- openFile f WriteMode   hq <- openFile q WriteMode@@ -108,16 +112,16 @@   hClose hq   hClose hf -hWriteFastaQual :: Handle -> Handle -> [Sequence] -> IO ()+hWriteFastaQual :: Handle -> Handle -> [Sequence a] -> IO () hWriteFastaQual hf hq = mapM_ wFQ     where wFQ s = wFasta hf s >> wQual hq s  -- | Lazily read sequence from handle-hReadFasta :: Handle -> IO [Sequence]+hReadFasta :: Handle -> IO [Sequence Unknown] hReadFasta h = (mkSeqs . B.lines) `fmap` B.hGetContents h  -- | Write sequences in FASTA format to a handle.-hWriteFasta :: Handle -> [Sequence] -> IO ()+hWriteFasta :: Handle -> [Sequence a] -> IO () hWriteFasta h = mapM_ (wFasta h)  wHead :: Handle -> SeqData -> IO ()@@ -126,17 +130,17 @@    B.hPut h l    B.hPut h $ B.pack "\n" -wFasta :: Handle -> Sequence -> IO ()+wFasta :: Handle -> Sequence a -> IO () wFasta h (Seq l d _) = do   wHead h l   let ls = splitsAt 60 d   mapM_ (B.hPut h) $ intersperse (B.pack "\n") ls   B.hPut h $ B.pack "\n" -hWriteQual :: Handle -> [Sequence] -> IO ()+hWriteQual :: Handle -> [Sequence a] -> IO () hWriteQual h = mapM_ (wQual h) -wQual :: Handle -> Sequence -> IO ()+wQual :: Handle -> Sequence a -> IO () wQual h (Seq l _ (Just q)) = do   wHead h l   let qls = splitsAt 20 q@@ -148,15 +152,15 @@ --   Blank lines are ignored. --   Comment lines start with "#" are allowed between sequences (and ignored). --   Lines starting with ">" initiate a new sequence.-mkSeqs :: [ByteString] -> [Sequence]+mkSeqs :: [ByteString] -> [Sequence Unknown] mkSeqs = map mkSeq . blocks -mkSeq :: [ByteString] -> Sequence+mkSeq :: [ByteString] -> Sequence Unknown mkSeq (l:ls) = Seq (B.drop 1 l) (B.concat $ takeWhile isSeq ls) Nothing     where isSeq s = (not . B.null) s && ((flip elem) (['A'..'Z']++['a'..'z']) . B.head) s mkSeq [] = error "empty input to mkSeq" -mkQual :: [ByteString] -> [Sequence]+mkQual :: [ByteString] -> [Sequence Unknown] mkQual = map f . blocks     where f (l:ls) = Seq (B.drop 1 l) B.empty                      (Just $ BB.pack (map int (B.words $ B.unlines ls)))

@@ -15,32 +15,33 @@ import System.IO  -- | Parse a .phd file, extracting the contents as a Sequence-readPhd :: FilePath -> IO Sequence+readPhd :: FilePath -> IO (Sequence Nuc) readPhd f = return . mkPhd =<< B.readFile f  -- | Parse .phd contents from a handle-hReadPhd :: Handle -> IO Sequence+hReadPhd :: Handle -> IO (Sequence Nuc) hReadPhd h = return . mkPhd =<< B.hGetContents h  -- | The actual phd parser.  --  aesthetics is not a major design goal... --  but error checking really should have been.  Sigh.-mkPhd :: B.ByteString -> Sequence+mkPhd :: B.ByteString -> (Sequence Nuc) mkPhd inp =    let (hd:fs) = filter (not . B.null) . B.lines $ inp       (comment,sd) = break (==B.pack "BEGIN_DNA") fs-      label = B.drop 15 hd+      (magic,label) = B.splitAt 15 hd+      more_magic = magic == B.pack "BEGIN_SEQUENCE "       fields = B.words . B.unlines                . filter (not . isSubstr (B.pack "_COMMENT")) $ comment       sdata = filter ((==3).length) . map B.words $ sd       err = error "failed to parse quality value"       qual = BB.fromChunks [BBB.pack . map (maybe err (fromIntegral . fst) . B.readInt . (!!1)) $ sdata]-  in qual `seq` (Seq (compact $ B.unwords (label:fields)) -      (compact $ B.concat $ map head sdata)-      (Just qual))---   Todo: check that we start with "BEGIN_SEQUENCE", and that we---   have a BEGIN_DNA/END_DNA region there.+  in if more_magic then qual `seq` (Seq (compact $ B.unwords (label:fields)) +                                    (compact $ B.concat $ map head sdata)+                                    (Just qual))+     else error "Incorrectly formatted PHD file - missing BEGIN_SEQUENCE"+     --   Todo: also check that we have a BEGIN_DNA/END_DNA region there.  isSubstr :: B.ByteString -> B.ByteString -> Bool isSubstr s = any (B.isPrefixOf s) . B.tails

@@ -13,16 +13,19 @@  module Bio.Sequence.SFF ( SFF(..), CommonHeader(..)                         , ReadHeader(..), ReadBlock(..)-                        , readSFF, writeSFF+                        , readSFF, writeSFF, writeSFF'                         , sffToSequence                         , test, convert-                        , Flow, Qual, Index+                        , Flow, Qual, Index, SeqData, QualData+                        , ReadName (..), decodeReadName, encodeReadName                         ) where  import Bio.Sequence.SeqData+import Bio.Sequence.SFF_name  import Data.Int import qualified Data.ByteString.Lazy as LB+import qualified Data.ByteString.Lazy.Char8 as LBC import qualified Data.ByteString as B import qualified Data.ByteString.Char8 as BC import Data.ByteString (ByteString)@@ -30,11 +33,12 @@  import Data.List (intersperse) import Data.Binary-import Data.Binary.Get (getByteString)+import Data.Binary.Get (getByteString,getLazyByteString) import qualified Data.Binary.Get as G-import Data.Binary.Put (putByteString)+import Data.Binary.Put (putByteString,putLazyByteString)  import Text.Printf (printf)+import System.IO  -- | The type of flowgram value type Flow = Int16@@ -51,17 +55,47 @@ readSFF :: FilePath -> IO SFF readSFF f = return . decode =<< LB.readFile f -sffToSequence :: SFF -> [Sequence]+sffToSequence :: SFF -> [Sequence Nuc] 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)++")!")+          clip (n,s,q) = let (k,s2) = LB.splitAt (fromIntegral $ key_length ch) s+                         in if k==lb (key ch) then Seq (lb n) s2 (Just $ LB.drop (fromIntegral $ key_length ch) q)+                                else error ("Couldn't match key in sequence "++BC.unpack n++" ("++LBC.unpack k++" vs. "++BC.unpack (key ch)++")!") +-- | Write an 'SFF' to the specified file name writeSFF :: FilePath -> SFF -> IO () writeSFF = encodeFile ++-- | Write an 'SFF' to the specified file name, but go back and+--   update the read count.  Useful if you want to output a lazy+--   stream of 'ReadBlock's.+writeSFF' :: FilePath -> SFF -> IO ()+writeSFF' f (SFF hs rs) = do+  h <- openFile f WriteMode+  LBC.hPut h $ encode hs+  c <- writeReads h (fromIntegral $ flow_length hs) rs+  putStrLn ("Count: "++show c)+  hSeek h AbsoluteSeek 20+  LBC.hPut h $ encode c+  hClose h++writeReads :: Handle -> Int -> [ReadBlock] -> IO Int32+writeReads _ _ [] = return 0+writeReads h i (r:rs) = do+  LBC.hPut h $ encode (RBI i r)+  c <- writeReads h i rs+  return (c+1)++data RBI = RBI Int ReadBlock++-- | Wrapper for ReadBlocks since they need additional information+instance Binary RBI where +    put (RBI c r) = do+      putRB c r+    get = undefined+       -- -------------------------------------------------- -- | test serialization by output'ing the header and first two reads  --   in an SFF, and the same after a decode + encode cycle.@@ -102,11 +136,12 @@                                    (do                                        rh <- get :: Get ReadHeader                                       let nb = fromIntegral $ num_bases rh+                                          nb' = fromIntegral $ num_bases rh                                           fl = fromIntegral $ flow_length chead                                       fg <- getByteString (2*fl)                                       fi <- getByteString nb-                                      bs <- getByteString nb-                                      qty <- getByteString nb+                                      bs <- getLazyByteString nb'+                                      qty <- getLazyByteString nb'                                       let l = (fl*2+nb*3) `mod` 8                                       when (l > 0) (skip (8-l))                                       return (ReadBlock rh (decodeArray fg) fi bs qty)@@ -115,7 +150,7 @@      put (SFF hd rds) = do       put hd-      mapM_ (putRB (fromIntegral $ flow_length hd)) rds+      mapM_ (put . RBI (fromIntegral $ flow_length hd)) rds  -- | A ReadBlock can't be an instance of Binary directly, since it depends on --   information from the CommonHeader.@@ -123,9 +158,11 @@ putRB fl rb = do   put (read_header rb)   mapM_ put (flowgram rb)+  -- ensure that flowgram has correct lenght+  replicateM (2*(fl-length (flowgram rb))) (put (0::Word8))   putByteString (flow_index rb)-  putByteString (bases rb)-  putByteString (quality rb)+  putLazyByteString (bases rb)+  putLazyByteString (quality rb)   let nb = fromIntegral $ num_bases $ read_header rb       l = (fl*2+nb*3) `mod` 8   when (l > 0) (pad (8-l))@@ -233,7 +270,9 @@       read_header                :: ReadHeader     -- The data block     , flowgram                   :: [Flow]-    , flow_index, bases, quality :: ByteString+    , flow_index                 :: ByteString+    , bases                      :: SeqData+    , quality                    :: QualData     }  instance Show ReadBlock where@@ -241,7 +280,7 @@         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+            , "bases:\t"++LBC.unpack b+            , "quality:\t"++(concat . intersperse " " . map show . LB.unpack) q             , ""             ])

@@ -0,0 +1,86 @@+module Bio.Sequence.SFF_name where++import qualified Data.ByteString.Char8 as B+import Data.ByteString.Char8 (ByteString, pack)+import Data.Array.Unboxed+import Data.Char (ord)++-- | Read names encode various information, as per this struct.+data ReadName = ReadName { date :: (Int,Int,Int)+                         , time :: (Int,Int,Int)+                         , region :: Int+                         , x_loc, y_loc :: Int } deriving Show++-- ----------------------------------------------------------+-- Decoding++decodeReadName :: ByteString -> Maybe ReadName+decodeReadName b = do t <- decodeDate $ B.take 6 b+                      r <- fst `fmap` (B.readInt $ B.take 2 $ B.drop 7 b)+                      l <- decodeLocation $ B.drop 9 b+                      return $ ReadName { date = (\[y,m,d] -> (y,m,d)) (take 3 t)+                               , time = (\[hh,mm,ss] -> (hh,mm,ss)) (drop 3 t)+                               , region = r+                               , x_loc = fst l, y_loc = snd l }++decodeLocation :: ByteString -> Maybe (Int,Int)+decodeLocation l = (`divMod` 4096) `fmap` decode36 l++decodeDate :: ByteString -> Maybe [Int]+decodeDate d    = (fixyear . reverse . (`divMods` [60,60,24,32,13])) =<< decode36 d+    where fixyear (i:is) = Just (2000+i:is)+          fixyear []     = Nothing++-- ----------------------------------------------------------+-- Encoding++encodeReadName :: ReadName -> ByteString+encodeReadName r =  B.concat [ encodeDate (date r) (time r) +                             , encodeRegion (region r)+                             , encodeLocation (x_loc r) (y_loc r)]++encodeLocation :: Int -> Int -> ByteString+encodeLocation = undefined++encodeRegion :: Int -> ByteString+encodeRegion = undefined++encodeDate :: (Int,Int,Int) -> (Int,Int,Int) -> ByteString+encodeDate = undefined++-- ----------------------------------------------------------++divMods :: Int -> [Int] -> [Int]+divMods x (i:is) = let (a,b) = x `divMod` i+                   in b : divMods a is+divMods x [] = [x]++-- ----------------------------------------------------------+-- Decoding base36 strings++decode36 :: ByteString -> Maybe Int+decode36 s = (foldr1 (\a b -> b*36+a) . reverse) `fmap` (mapM decCh . B.unpack $ s)++{-+decode36' = dec 0+    where dec i b = case uncons b of Just (c,rest) -> dec (i*36+fromJust (decCh c)) rest+                                     Nothing       -> i+          fromJust (Just z) = z+-}++decCh :: Char -> Maybe Int+decCh x | x >= 'A' && x <= 'Z' = Just (ord x - ord 'A')+        | x >= '0' && x <= '9' = Just (26 + ord x - ord '0')+        | otherwise            = Nothing -- error ("decode36: can't decode "++show x)++encode36 :: Int -> ByteString+encode36 = pack . map (b36!) . reverse . enc+    where+      enc 0 = []+      enc i = let (a,b) = i `divMod` 36+              in b : enc a++b36 :: UArray Int Char+b36 = listArray (0,35) (['A'..'Z']++['0'..'9'])++

@@ -5,10 +5,14 @@    like @revcompl@, only makes sense for nucleotides. -} +{-# LANGUAGE EmptyDataDecls #-}+ module Bio.Sequence.SeqData      (       -- * Data structure-      -- | A sequence is a header, sequence data itself, and optional quality data.+      -- | A sequence is a header, sequence data itself, and optional quality data. +      --   Sequences are type-tagged to identify them as nucleotide, amino acids,+      --   or unknown type.        --   All items are lazy bytestrings.  The Offset type can be used for indexing.       Sequence(..), Offset, SeqData, @@ -25,25 +29,37 @@       -- * Converting to and from [Char]     , fromStr, toStr +      -- * Sequence utilities+    , defragSeq+       -- * Nucleotide functionality       -- | Nucleotide sequences contain the alphabet [A,C,G,T].       -- IUPAC specifies an extended nucleotide alphabet with wildcards, but       -- it is not supported at this point.-    , compl, revcompl+    , compl, revcompl, revcompl', Nuc, castToNuc        -- * Protein functionality       -- | Proteins are chains of amino acids, represented by the IUPAC alphabet.     , Amino(..), translate, fromIUPAC, toIUPAC -- amino acids+    , castToAmino+    +      -- * Display a nicely formated sequence.+    , putSeqLn, seqToStr++      -- * Default type for sequences+    , Unknown     ) where -import Data.List (unfoldr)-import Data.Char (toUpper)+import Data.List (unfoldr, intercalate, isPrefixOf)+import Data.Char (toUpper, isNumber) import Data.Int import Data.Word import Data.Maybe (fromJust, isJust) import qualified Data.ByteString.Lazy.Char8 as B import qualified Data.ByteString.Lazy as BB-+import qualified Data.ByteString.Char8 as BS+import qualified Data.ByteString as BBS+import Text.Printf (printf)  -- | An offset, index, or length of a 'SeqData' type Offset  = Int64@@ -59,9 +75,186 @@ type QualData = BB.ByteString -- mild abuse  -- | A sequence consists of a header, the sequence data itself, and optional quality data.-data Sequence = Seq !SeqData !SeqData !(Maybe QualData) -- ^ header and actual sequence-                deriving (Show,Eq)+--   The type parameter is a phantom type to separate nucleotide and amino acid sequences+data Sequence t = Seq !SeqData !SeqData !(Maybe QualData) -- ^ header and actual sequence+                  deriving Eq +-- | For type tagging sequences (protein sequences use 'Amino' below)+data Nuc+data Unknown++-- | Phantom type functionality+castToNuc :: Sequence a -> Sequence Nuc+castToNuc (Seq h d mq) = Seq h d mq++castToAmino :: Sequence a -> Sequence Amino+castToAmino (Seq h d mq) = Seq h d mq+++-- | A more arranged show instance for Sequences reassembling the display of+--   the fasta-format+instance Show (Sequence a) where+    show s = seqToStr s 6 10 []++-- | Returns a properly formatted and probably highlighted string+-- | representation of a sequence. Highlighting is done using ANSI-Escape+-- | sequences.+seqToStr :: Sequence a -> Int -> Int -> [(Int, Int)] -> [Char]+seqToStr (Seq header sd _) chunks len parts =+    let strSeq     = toStr sd+        strHeader  = toStr header+        identifier = head (words strHeader)+        comment    = intercalate " " $ tail (words strHeader)++        numLength  = (length . show . length) strSeq+        width      = chunks * len + numLength + chunks+        format     = intercalate "\n" . splits (width)+    in+    (line "ID" width)      ++ (format identifier) ++ "\n\n" +++    (line "COMMENT" width) ++ (format comment)    ++ "\n\n" +++    (line "DATA" width)    ++ (showDNA strSeq len chunks parts)+ +-- | A simple function to display a sequence: we generate the sequence string and+-- | call putStrLn+putSeqLn :: Sequence a -> Int -> Int -> [(Int, Int)] -> IO ()+putSeqLn s chunks len parts = putStrLn $ seqToStr s chunks len parts+ +-- Creates a header line of length n+line :: [Char] -> Int -> [Char]+line header width =+  header ++ " " ++ (take (width - (length header) - 1) (repeat '-')) ++ "\n"+ +-- | Splits a string into parts of size width. The last element can be shorter.+splits :: Int -> [t] -> [[t]]+splits _     [] = []+splits width xs =+    let (a, b) = splitAt width xs in+    [a] ++ (splits width b)++-- Highlighting and Dehighlighting ANSI-Escape sequences+beginHighlight, endHighlight :: String+beginHighlight = "\ESC[7m"+endHighlight   = "\ESC[0m"++-- Creates a string representation of the dna-part of a sequence. We can either+-- highlight specific areas, or, if this is empty, numberize the rows. At the+-- moment, both at the same time is not possible (see comment in the numberize-+-- function). If you want numberized output, take a look at the showDNANumber-+-- function.+showDNA :: [Char] -> Int -> Int -> [(Int, Int)] -> [Char]+showDNA xs len chunks highs =+    let hs  = highlight xs highs+        ks  = splitsSkip hs len [beginHighlight, endHighlight]+        ch  = chunkify ks chunks       +        ns  = numberize ch chunks len+        ns' = clean ns False+    in ns'+    +-- Add orientation numbers to each line.+numberize :: [[Char]] -> Int -> Int -> [Char]+numberize ls chunks len = +    let -- The initial numbering of the single elements in lines+        numbered = zip [0, (chunks * len)..] ls +        +        -- Create the format string with the right alignment, depending on the +        -- number of elements in lines+        numLength   = (show . length . show ) (chunks * len * (length ls))+        formatStr   = "%" ++ numLength ++ "d "+        toNum (a,b) = (printf formatStr a) ++ b        +    in intercalate "\n" (map toNum numbered)+  +-- Concatenates each n elements of a list.    +chunkify :: [[Char]] -> Int -> [[Char]]+chunkify [] _n = []+chunkify xs  n =  +    [intercalate " " (take n xs)] ++ (chunkify (drop n xs) n)+      +-- Highlighting should not occur on spaces or over newline borders. We walk+-- through a generated string and add appropiate stop-highlight and start-+-- highlight markers.+--+-- This is one of the functions which can be definitely optimized by a seasoned+-- Haskell programmer.+clean :: [Char] -> Bool -> [Char]+clean []         _      = []++-- Handling spaces+clean (' ':xs)   True   = endHighlight ++ " " ++ beginHighlight ++ clean xs True+clean (' ':xs)   False  = " " ++ clean xs False++-- Newlines and numbering+clean ('\n':xs)   True   = +    let (num, dna) = span (\c -> isNumber c) xs+    in endHighlight ++ "\n" ++ num ++ beginHighlight ++ clean dna True+clean ('\n':xs)   False  = "\n" ++ clean xs False++-- Checking for highlighting at all+clean str@(x:xs) isHigh =+    -- Special cases: check for the beginning or ending of highlighting +    let a = isPrefixOf beginHighlight str+        b = isPrefixOf endHighlight str+    in if a then x : clean xs True+            else if b then x : clean xs False+                      else x : clean xs isHigh++-- Highlights a range of strings. Each start and end to be highlighted is +-- described by a position-tuple (start, end).+highlight :: [Char] -> [(Int, Int)] -> [Char]+highlight ss [] = ss +highlight ss ((start, end):rest) = +    let sl = [beginHighlight, endHighlight]+        a  = insert ss start   beginHighlight sl+        b  = insert a  (end+1) endHighlight   sl+    in highlight b rest++-- Generates lists of length len from xs, skipping sublist in skips while+-- counting list elements.+splitsSkip :: (Eq t) => [t] -> Int -> [[t]] -> [[t]]+splitsSkip [] _   _     = []+splitsSkip xs len skips =+    (takeSkip xs len skips) : splitsSkip (dropSkip xs len skips) len skips+    +-- Emulates the normal take behaviour, skipping over sublists in skips.+takeSkip :: (Eq a) => [a] -> Int -> [[a]] -> [a]+takeSkip []         _ _     = []+takeSkip _          0 _     = []+takeSkip str@(x:xs) n skips =+    let n' = case (testPrefixes str skips) of+               Nothing  -> n-1+               Just len -> n + len - 1+    in x : takeSkip xs n' skips++-- Emulates the normal drop behaviour, skipping over sublists in skips.+dropSkip :: (Eq a) => [a] -> Int -> [[a]] -> [a]+dropSkip []         _ _     = []+dropSkip xs         0 _     = xs+dropSkip str@(x:xs) n skips =+    let n' = case (testPrefixes str skips) of+               Nothing  -> n-1+               Just len -> n + len - 1+    in dropSkip xs n' skips++-- Another definition of insert, which provides support for handling skipping+-- of strings. Inserts elem *before* pos into xs, but skips all strings in+-- skips while calculating the specified position+insert :: (Eq a) => [a] -> Int -> [a] -> [[a]] -> [a]+insert []         _   _    _     = []+insert xs         0   elem _     = elem ++ xs+insert str@(x:xs) pos elem skips = +    let pos' = case (testPrefixes str skips) of+               Nothing  -> pos-1+               Just len -> pos + len - 1+    in x : insert xs pos' elem skips+        +-- Test, if one of the prefixes is the prefix of xs. If yes, return it's length+-- else Nothing  +testPrefixes :: (Eq a) => [a] -> [[a]] -> Maybe Int+testPrefixes _  []     = Nothing+testPrefixes xs (p:ps) = +    if (isPrefixOf p xs)+        then return (length p)+        else testPrefixes xs ps+         -- | Convert a String to 'SeqData' fromStr :: String -> SeqData fromStr = B.pack@@ -72,47 +265,63 @@  -- | Read the character at the specified position in the sequence. {-# INLINE (!) #-}-(!) :: Sequence -> Offset -> Char+(!) :: Sequence a -> Offset -> Char (!) (Seq _ bs _) = B.index bs  {-# INLINE (?) #-}-(?) :: Sequence -> Offset -> Qual+(?) :: Sequence a -> Offset -> Qual (?) (Seq _ _ (Just qs)) = BB.index qs (?) (Seq _ _ Nothing) =      error "Attempting to use 'seqqual' on sequence without associated quality data"  -- | Return sequence length.-seqlength :: Sequence -> Offset+seqlength :: Sequence a -> Offset seqlength (Seq _ bs _) = B.length bs  -- | Return sequence label (first word of header)-seqlabel :: Sequence -> SeqData+seqlabel :: Sequence a -> SeqData seqlabel (Seq l _ _) = head (B.words l)  -- | Return full header.-seqheader :: Sequence -> SeqData+seqheader :: Sequence a -> SeqData seqheader (Seq l _ _) = l  -- | Return the sequence data.-seqdata :: Sequence -> SeqData+seqdata :: Sequence a -> SeqData seqdata (Seq _ d _) = d  -- | Return the quality data, or error if none exist.  Use hasqual if in doubt.-seqqual :: Sequence -> QualData+seqqual :: Sequence a -> QualData seqqual (Seq _ _ (Just q)) = q seqqual (Seq _ _ Nothing) =      error "Attempting to use 'seqqual' on sequence without associated quality data"  -- | Check whether the sequence has associated quality data.-hasqual :: Sequence -> Bool+hasqual :: Sequence a -> Bool hasqual (Seq _ _ q) = isJust q  -- | Modify the header by appending text, or by replacing --   all but the sequence label (i.e. first word).-appendHeader, setHeader :: Sequence -> String -> Sequence+appendHeader, setHeader :: Sequence a -> String -> Sequence a appendHeader (Seq h d q) t = (Seq (B.append h (B.pack (" "++t))) d q) setHeader (Seq h d q) t = (Seq (B.append (head $ B.words h) (B.pack (" "++t))) d q) +-- | Returns a sequence with all internal storage freshly copied and+-- with sequence and quality data present as a single chunk.  +--  +-- By freshly copying internal storage, 'defragSeq' allows garbage+-- collection of the original data source whence the sequence was+-- read; otherwise, use of just a short sequence name can cause an+-- entire sequence file buffer to be retained.+-- +-- By compacting sequence data into a single chunk, 'defragSeq' avoids+-- linear-time traversal of sequence chunks during random access into+-- sequence data.+defragSeq :: Sequence t -> Sequence t+defragSeq (Seq name sequ qual) = Seq (defragB name) (defragB sequ) (fmap defragBB qual)+    where defragB = B.fromChunks . (: []) . BS.concat . B.toChunks+          defragBB = BB.fromChunks . (: []) . BBS.concat . BB.toChunks+ ------------------------------------------------------------ -- Nucleotide (DNA, RNA) specific stuff ------------------------------------------------------------@@ -120,10 +329,14 @@ -- | Calculate the reverse complement. --   This is only relevant for the nucleotide alphabet,  --   and it leaves other characters unmodified.  -revcompl :: Sequence -> Sequence-revcompl (Seq l bs qs) = Seq l (B.reverse $ B.map compl bs)+revcompl :: Sequence Nuc -> Sequence Nuc+revcompl (Seq l bs qs) = Seq l (revcompl' bs)                         (maybe Nothing (Just . BB.reverse) qs) +-- | Calculate the reverse complent for SeqData only.+revcompl' :: SeqData -> SeqData+revcompl' = B.reverse . B.map compl+ -- | Complement a single character.  I.e. identify the nucleotide it  --   can hybridize with.  Note that for multiple nucleotides, you usually --   want the reverse complement (see 'revcompl' for that).@@ -145,7 +358,7 @@ -- | Translate a nucleotide sequence into the corresponding protein --   sequence.  This works rather blindly, with no attempt to identify ORFs --   or otherwise QA the result.-translate :: Sequence -> Offset -> [Amino]+translate :: Sequence Nuc -> Offset -> [Amino] translate s' o' = unfoldr triples (s',o')    where triples (s,o) =               if o > seqlength s - 3 then Nothing

@@ -1,4 +1,5 @@ {- |+    This module implements the 2bit format for sequences.     Based on: <http://genome.ucsc.edu/FAQ/FAQformat#format7>@@ -11,9 +12,15 @@   module Bio.Sequence.TwoBit-   ( decode2Bit, read2Bit, hRead2Bit+   ( decode2Bit, +     read2Bit, +     hRead2Bit,+     encode2Bit,+     write2Bit,+     hWrite2Bit    ) where + import Bio.Sequence.SeqData import qualified Data.ByteString.Lazy as BB import qualified Data.ByteString.Lazy.Char8 as B@@ -30,13 +37,13 @@ import Test.QuickCheck hiding (check)    -- QC 1.0 -- import Test.QuickCheck hiding ((.&.)) -- QC 2.0 -type ByteString = B.ByteString  -- constants default_magic, default_version :: Word32 default_magic   = 0x1A412743 default_version = 0 + -- binary extras check :: Monad m => (a -> Bool) ->  a -> m a check p x = if (p x) then return x else fail "check failed"@@ -51,6 +58,8 @@ unbytes :: Integral a => [Word8] -> a unbytes = Data.List.foldr (\x y -> y*256+x) 0 . map fromIntegral ++ -- Conflicts with Bio.Util.TestBase -- instance Arbitrary Word8 where --    arbitrary = choose (0,255::Int) >>= return . fromIntegral@@ -58,166 +67,391 @@ -- 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+++++-- "in-core" representation of 2Bit data types+ data Header = Header { swap :: Bool,  version, count, reserved :: Word32 } instance Show Header where show (Header _ v c r) = "H "++show (v,c,r)  instance Binary Header where     get = do-      m <- get-      v <- get >>= check (==default_version)-      c <- get-      r <- get-      let s = if m==default_magic then id else if m==bswap 4 default_magic then bswap 4-                                               else error "2bit decode: incorrect magic number"-      return (Header (m/=default_magic) (s v) (s c) (s r))+       m <- get+       v <- get >>= check (==default_version)+       c <- get+       r <- get+       let s = if m == default_magic +                   then id +                   else if m == bswap 4 default_magic +                            then bswap 4+                            else error "2bit decode: incorrect magic number"+       return (Header (m /= default_magic) (s v) (s c) (s r))     put (Header m v c r) = do-      put default_magic-      put default_version-      put c-      put (0 :: Word32)+       put default_magic+       put default_version+       put c+       put (0 :: Word32) -data Entry  = Entry  { name :: ByteString, offset :: Word32 }+++data Entry  = Entry  { name :: B.ByteString, offset :: Word32 }               deriving Show -swapEntry e = e { offset = bswap 4 (offset e) }+-- Byte swap an Entry's offset!+swapEntry :: Entry -> Entry+swapEntry entry = entry { offset = bswap 4 (offset entry) }  instance Binary Entry where     get = do-      len <- getWord8-      name <- replicateM (fromIntegral len) getWord8-      offset <- get-      return (Entry (BB.pack name) offset)-    put (Entry bs offset) = do-      let l = fromIntegral $ B.length bs :: Word8-      put l-      mapM_ put $ BB.unpack bs-      put offset+       len <- getWord8+       name <- replicateM (fromIntegral len) getWord8+       offset <- get+       return (Entry (BB.pack name) offset)+    put (Entry byteString offset) = do+       let len = fromIntegral $ B.length byteString :: Word8+       put len  +       mapM_ put $ BB.unpack byteString+       put offset ++ data Entries = Entries Header [Entry]-instance Show Entries where show (Entries h es) = unlines (show h : map show es) +instance Show Entries +   where show (Entries h es) = unlines (show h : map show es)+ instance Binary Entries where     get = do-      h <- get-      es <- replicateM (fromIntegral $ count h) get-      return (Entries h $ if swap h then map swapEntry es else es)+       h <- get+       es <- replicateM (fromIntegral $ count h) get+       return (Entries h $ if swap h then map swapEntry es else es)+    put (Entries h es) = do+       put h+       put es    ++ {-- dnaSize - number of bases of DNA in the sequence- nBlockCount - the number of blocks of Ns in the file (representing unknown sequence)- nBlockStarts - the starting position for each block of Ns- nBlockSizes - the size of each block of Ns- maskBlockCount - the number of masked (lower-case) blocks++   Sequence Record definition++ dnaSize         - number of bases of DNA in the sequence+ nBlockCount     - the number of blocks of Ns in the file (representing unknown sequence)+ nBlockStarts    - the starting position for each block of Ns+ nBlockSizes     - the size of each block of Ns+ maskBlockCount  - the number of masked (lower-case) blocks  maskBlockStarts - the starting position for each masked block- maskBlockSizes - the size of each masked block- packedDna - the DNA packed to two bits per base+ maskBlockSizes  - the size of each masked block+ packedDna       - the DNA packed to two bits per base+ -} -data SR = SR { dnaSize, nBlockCount :: Word32-                         , nBlockStarts, nBlockSizes :: [Word32]-                         , maskBlockCount :: Word32-                         , maskBlockStarts, maskBlockSizes :: [Word32]-                         , packedDna :: [Word8]-                         , reserved2 :: Word32 }++data SR = SR { dnaSize         :: Word32,+               nBlockCount     :: Word32,+               nBlockStarts    :: [Word32],+               nBlockSizes     :: [Word32],+               maskBlockCount  :: Word32,+               maskBlockStarts :: [Word32],+               maskBlockSizes  :: [Word32],+               packedDna       :: [Word8],+               reserved2       :: Word32 }         deriving Show ++ -- big- and little-endian variants (what a mess) newtype SRBE = SRBE SR deriving Show newtype SRLE = SRLE SR deriving Show + instance Binary SRBE where-    get = undefined+    get = do+       dz <- get :: Get Word32+       nc <- get :: Get Word32+       let n = fromIntegral nc+       nbs <- replicateM n get+       nbsz <- replicateM n get+       mc <- get :: Get Word32+       let m = fromIntegral mc+       mbs <- replicateM m get+       mbsz <- replicateM m get+       _reserved <- get :: Get Word32 -- !!!! oops?+       let d = fromIntegral dz+       pdna <- replicateM ((d+3) `div` 4)  get+       return (SRBE $ SR dz nc nbs nbsz+              mc mbs mbsz pdna _reserved)+    -- should this happen?  Why not just write default format?+    put (SRBE sr) = do+       put $ dnaSize sr+       put $ nBlockCount sr+       mapM_ put (nBlockStarts sr)+       mapM_ put (nBlockSizes sr)+       put $ maskBlockCount sr+       mapM_ put (maskBlockStarts sr)+       mapM_ put (maskBlockSizes sr)+       put (0::Word32)+       mapM_ put (packedDna sr) + instance Binary SRLE where     get = do-      dz <- get :: Get Word32-      nc <- get :: Get Word32-      let n = fromIntegral $ bswap 4 nc-      nbs <- replicateM n get-      nbsz <- replicateM n get-      mc <- get :: Get Word32-      let m = fromIntegral $ bswap 4 mc-      mbs <- replicateM m get-      mbsz <- replicateM m get-      _reserved <- get :: Get Word32 -- !!!! oops?-      let d = fromIntegral $ bswap 4 dz-      pdna <- replicateM ((d+3) `div` 4)  get-      return (SRLE $ SR (bswap 4 dz) (bswap 4 nc)+       dz <- get :: Get Word32+       nc <- get :: Get Word32+       let n = fromIntegral $ bswap 4 nc+       nbs <- replicateM n get+       nbsz <- replicateM n get+       mc <- get :: Get Word32+       let m = fromIntegral $ bswap 4 mc+       mbs <- replicateM m get+       mbsz <- replicateM m get+       _reserved <- get :: Get Word32 -- !!!! oops?+       let d = fromIntegral $ bswap 4 dz+       pdna <- replicateM ((d+3) `div` 4)  get+       return (SRLE $ SR (bswap 4 dz) (bswap 4 nc)               (map (bswap 4) nbs) (map (bswap 4) nbsz)               (bswap 4 mc) (map (bswap 4) mbs) (map (bswap 4) mbsz) pdna (bswap 4 _reserved))     -- should this happen?  Why not just write default format?     put (SRLE sr) = do-      put (bswap 4 $ dnaSize sr)-      put (bswap 4 $ nBlockCount sr)-      mapM_ (put . bswap 4) (nBlockStarts sr)-      mapM_ (put . bswap 4) (nBlockSizes sr)-      put (bswap 4 $ maskBlockCount sr)-      mapM_ (put . bswap 4) (maskBlockStarts sr)-      mapM_ (put . bswap 4) (maskBlockSizes sr)-      put (0::Word32)-      mapM_ put (packedDna sr)+       put (bswap 4 $ dnaSize sr)+       put (bswap 4 $ nBlockCount sr)+       mapM_ (put . bswap 4) (nBlockStarts sr)+       mapM_ (put . bswap 4) (nBlockSizes sr)+       put (bswap 4 $ maskBlockCount sr)+       mapM_ (put . bswap 4) (maskBlockStarts sr)+       mapM_ (put . bswap 4) (maskBlockSizes sr)+       put (0::Word32)+       mapM_ put (packedDna sr) --- used to convert to/from sequence data in the Sequence data structure-fromSR :: SR -> ByteString+++-- Used to convert from sequence data in the Sequence data structure to ByteString+fromSR :: SR -> B.ByteString fromSR sr = B.unfoldr go (0,low,ns,take (fromIntegral $ dnaSize sr) dna)     where-    low = combine (maskBlockStarts sr) (maskBlockSizes sr)-    ns  = combine (nBlockStarts sr) (nBlockSizes sr)-    combine starts lengths = concatMap (\(p,l) -> [p..p+l-1]) $ zip starts lengths+       low = combine (maskBlockStarts sr) (maskBlockSizes sr)+       ns  = combine (nBlockStarts sr) (nBlockSizes sr) -    dna = decodeDNA $ packedDna sr-    decodeDNA = concatMap (\x -> [shiftR (x .&. 0xC0) 6, shiftR (x .&. 0x30) 4, shiftR (x .&. 0x0C) 2, x .&. 0x03])-    dec1 x = case x of 0 -> 'T'; 1 -> 'C'; 2 -> 'A'; 3 -> 'G'; _ -> error ("can't decode value "++show x)+       combine :: (Num t, Enum t) => [t] -> [t] -> [t]+       combine starts lengths = concatMap (\(p,l) -> [p..p+l-1]) $ zip starts lengths -    go (_,_,_,[]) = Nothing-    go (pos,(l:ls),(n:ns),(d:ds))-        | pos == l && pos == n = Just ('n',(pos+1,ls,ns,ds))-        | pos == l             = Just (toLower (dec1 d),(pos+1,ls,n:ns,ds))-        |             pos == n = Just ('N',(pos+1,l:ls,ns,ds))-        | otherwise            = Just (dec1 d, (pos+1,l:ls,n:ns,ds))-    go (pos,[],n:ns,d:ds)-        |             pos == n = Just ('N',(pos+1,[],ns,ds))-        | otherwise            = Just (dec1 d, (pos+1,[],n:ns,ds))-    go (pos,l:ls,[],d:ds)-        | pos == l             = Just (toLower (dec1 d),(pos+1,ls,[],ds))-        | otherwise            = Just (dec1 d, (pos+1,l:ls,[],ds))-    go (pos,[],[],d:ds)        = Just (dec1 d, (pos+1,[],[],ds))+       -- Unpack a 2Bit packed DNA sequence into a list of 2 bit values 0-3 +       dna = decodeDNA $ packedDna sr+       decodeDNA = concatMap (\x -> [shiftR (x .&. 0xC0) 6, shiftR (x .&. 0x30) 4, shiftR (x .&. 0x0C) 2, x .&. 0x03])+++       -- Map 2Bit nucleotide encodings to their character equivalents+       dec1 :: (Num t) => t -> Char+       dec1 x = case x of +                   0 -> 'T'; +                   1 -> 'C'; +                   2 -> 'A'; +                   3 -> 'G'; +                   _ -> error ("can't decode value "++show x)+++       go :: (Num a, Num t) => (a, [a], [a], [t]) -> Maybe (Char, (a, [a], [a], [t]))+       go (_,_,_,[])              = Nothing+       go (pos,(l:ls),(n:ns),(d:ds))+           | pos == l && pos == n = Just ('n',(pos+1,ls,ns,ds))+           | pos == l             = Just (toLower (dec1 d),(pos+1,ls,n:ns,ds))+           |             pos == n = Just ('N',(pos+1,l:ls,ns,ds))+           | otherwise            = Just (dec1 d, (pos+1,l:ls,n:ns,ds))+       go (pos,[],n:ns,d:ds)+           |             pos == n = Just ('N',(pos+1,[],ns,ds))+           | otherwise            = Just (dec1 d, (pos+1,[],n:ns,ds))+       go (pos,l:ls,[],d:ds)+           | pos == l             = Just (toLower (dec1 d),(pos+1,ls,[],ds))+           | otherwise            = Just (dec1 d, (pos+1,l:ls,[],ds))+       go (pos,[],[],d:ds)        = Just (dec1 d, (pos+1,[],[],ds)) --    go x = error (show x) +++toSR :: B.ByteString -> SR+toSR bs = undefined++++splits :: [Int64] -> B.ByteString -> [B.ByteString] +splits [] cs = [cs]+splits (e:es) cs = let (this,rest) = B.splitAt e cs+                   in this : splits es rest++++ -- | Parse a (lazy) ByteString as sequences in the 2bit format.-decode2Bit :: B.ByteString -> [Sequence]-decode2Bit cs = let (Entries h es) = decode cs :: Entries+decode2Bit :: B.ByteString -> [Sequence Unknown]+decode2Bit cs = let +                    -- decode to (Header, [Entry]) from ByteString+                    (Entries h es) = decode cs :: Entries++                    -- map to [64-bit offset] from [32-bit offset]!!!                      ms = map (fromIntegral . offset) es++                    -- break the ByteString up into chunks of raw 2Bit "sequences"                     (c:chunks) = zipWith (-) ms (0:ms)++                    --  build raw undecoded [Sequence]!!!+                    --  chunks :: [Int64]+                    --          "drop c cs" .. chop off all cruft in the beginning!!!                     ss = splits chunks $ B.drop c cs-               in map ($ Nothing) $ zipWith Seq (map name es)-                      $ map fromSR $ case swap h of-                                       True -> map (unSRLE.decode) ss-                                       False -> map (unSRBE.decode) ss -splits [] cs = [cs]-splits (e:es) cs = let (this,rest) = B.splitAt e cs-                   in this : splits es rest +                --                              ($ Nothing) :: (Maybe a -> b) -> b+                --                              zipWith :: (a -> b -> c) -> [a] -> [b] -> [c]+                --                              Seq :: SeqData -> SeqData -> Maybe QualData -> Sequence+                --                              map fromSR :: [SR] -> [B.ByteString] +                --                              map (unSRLE.decode) :: [B.ByteString] -> [SR]+                --                              map (unSRBE.decode) :: [B.ByteString] -> [SR]+                in map ($ Nothing) $ zipWith Seq (map name es) $ map fromSR $ case swap h of+                                                                                 True -> map (unSRLE.decode) ss+                                                                                 False -> map (unSRBE.decode) ss+                   ++++type SequenceLabel = SeqData+type SequenceSize = Offset+type SequenceData = SeqData +type TwoBitData = (SequenceLabel, SequenceSize, SequenceData)++-- | Marshall from neutral representation to the 2Bit ByteString rep+encode2Bit :: [Sequence a] -> B.ByteString+encode2Bit ss = let++                   buildHeader :: [Sequence a] -> Header+                   buildHeader ss = Header {swap = True, +                                            version = default_version, +                                            count = sequenceListLength ss, +                                            reserved = 0}+++                   -- Build the list of 2Bit Entries +                   buildEntries :: [TwoBitData] -> Offset -> [Entry]+                   buildEntries [] _ = []+                   buildEntries ((label, length, _):xs) currentOffset = Entry {name=label, +                                                                               offset=fromIntegral $ currentOffset} : buildEntries xs (currentOffset+length)   ++                   -- Build a 2Bit Sequence Record!!+                   buildSR :: TwoBitData -> SR+                   buildSR (_, size, dnaData) = SR {dnaSize = (fromIntegral size), +                                                    nBlockCount = 0,+                                                    nBlockStarts = [],+                                                    nBlockSizes = [],+                                                    maskBlockCount = 0,+                                                    maskBlockStarts = [], +                                                    maskBlockSizes = [],+                                                    packedDna = encodeDNA $ splitWord8 $ explode dnaData, +                                                    reserved2 = 0}+++                   -- Total # of sequences present+                   sequenceListLength :: [Sequence a] -> Word32+                   sequenceListLength [] = 0+                   sequenceListLength (s:ss) = 1 + sequenceListLength ss+  ++                   -- Build a list of vital data to do 2Bit marshalling+                   sequenceListExtract :: [Sequence a] -> [TwoBitData]+                   sequenceListExtract ss = map (\seq -> (seqlabel seq, seqlength seq, seqdata seq)) ss++  +                   -- Map a nucleotide to its respective 2Bit encoding+                   enc1 :: (Num t) => Char -> t +                   enc1 c = case c of +                               'T' -> 0; +                               'C' -> 1; +                               'A' -> 2; +                               'G' -> 3;+++                   -- Take a ByteString of nucleotides in character encoding and map to a list of their corresponding TwoBit encodings+                   explode :: SeqData -> [Word8]+                   explode seq = map enc1 $ B.unpack seq +++                   -- Split a list of Word8's into 4 character sublists (i.e. quads) +                   splitWord8 ::  [Word8] -> [[Word8]]+                   splitWord8 [] = []+                   splitWord8 cs = let (this,rest) = splitAt 4 cs+                                      in this : splitWord8 rest+++                   -- Build a 2Bit packed DNA sequence+                   encodeDNA :: [[Word8]] -> [Word8]+                   encodeDNA [] = []+                   encodeDNA (w:ws) = (pack2Bit $ getQuad) : encodeDNA ws++                                     where++                                        -- Length of a sequence+                                        len = length w+++                                        -- Build a 2Bit encoded Word8 from "exploded" encoding +                                        pack2Bit :: (Word8, Word8, Word8, Word8) -> Word8+                                        pack2Bit (d1, d2, d3, d4) = shiftL (d1) 6 .|. shiftL (d2) 4 .|. shiftL (d3) 2 .|. d4++                                        -- Build a quad, i.e. a 4 tuple of Word8+                                        getQuad :: (Word8, Word8, Word8, Word8)+                                        getQuad = case len of+                                                        1 -> (w !! 0, 0, 0, 0);+                                                        2 -> (w !! 0, w !! 1, 0, 0);+                                                        3 -> (w !! 0, w !! 1, w !! 2, 0);+                                                        4 -> (w !! 0, w !! 1, w !! 2, w !! 3);++                          ++               in ++                   -- Serialize/marshall into ByteString representation+                   B.append (encode (buildHeader ss))+                      (B.append+                         (B.concat (map encode (buildEntries (sequenceListExtract ss) 56)))  -- TEMP WNH+                         (B.concat (map (encode . SRBE) (map buildSR (sequenceListExtract ss)))))      -- TEMP WNH+                         -- WNH (B.concat (map (encode . SRBE) (map buildSR (sequenceListExtract ss)))))++++unSRBE :: SRBE -> SR unSRBE (SRBE x) = x++unSRLE :: SRLE -> SR unSRLE (SRLE x) = x -toSR :: ByteString -> SR-toSR bs = undefined --- | Extract sequences from a file in 2bit format.-read2Bit  :: FilePath -> IO [Sequence]+++-- | Read sequences from a file in 2bit format and +-- | unmarshall/deserialize into Sequence format.+read2Bit  :: FilePath -> IO [Sequence Unknown] read2Bit f = B.readFile f >>= return . decode2Bit --- | Extract sequences in the 2bit format from a handle.-hRead2Bit :: Handle   -> IO [Sequence]+-- | Read sequences from a file handle in the 2bit format and+-- | unmarshall/deserialze into Sequence format.+hRead2Bit :: Handle   -> IO [Sequence Unknown] hRead2Bit h = B.hGetContents h >>= return . decode2Bit --- | Write sequences to file in the 2bit format.-write2Bit  :: FilePath -> [Sequence] -> IO ()-write2Bit = undefined --- | Write sequences to a handle in the 2bit format.-hWrite2Bit :: Handle   -> [Sequence] -> IO ()-hWrite2Bit = undefined+++-- | Marshall/serialize [Sequence] into 2Bit format and write to a file. +write2Bit  :: FilePath -> [Sequence a] -> IO ()+write2Bit f seq = do+                   let byteString = encode2Bit seq    +                   B.writeFile f byteString+                   return ()++-- | Marshall/serialize [Sequence] into 2Bit format and write to a file using handle. +hWrite2Bit :: Handle   -> [Sequence a] -> IO ()+hWrite2Bit h seq = do+                   let byteString = encode2Bit seq+                   B.hPut h byteString+                   return ()+++

@@ -25,14 +25,14 @@ 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+newtype EST = E (Sequence Nuc) deriving Show+newtype ESTq = Eq (Sequence Nuc) deriving Show+newtype Protein = P (Sequence Amino) 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+newtype EST_short = ES (Sequence Nuc) deriving Show+newtype EST_long  = EL (Sequence Nuc) deriving Show+newtype EST_set  = ESet [Sequence Nuc] deriving Show  -- | Take time (CPU and wall clock) and report it time :: String -> IO () -> IO ()

@@ -1,5 +1,5 @@ Name:                bio-Version:             0.3.5+Version:             0.4 License:             LGPL License-file:        LICENSE Author:              Ketil Malde@@ -8,15 +8,16 @@ Category:            Bioinformatics Synopsis:            A bioinformatics library Description:         This is a collection of data structures and algorithms-                     I've found useful when building various bioinformatics-related tools+                     useful for building bioinformatics-related tools                      and utilities.                      .                      Current list of features includes: a Sequence data type supporting-                     protein and nucleotide sequences and conversion between them, quality+                     protein and nucleotide sequences and conversion between them.  As of version+		     0.4, different kinds of sequence have different types.  Support for quality                      data, reading and writing Fasta formatted files, reading TwoBit and-                     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+                     phd formats, and Roche/454 SFF files.  Rudimentary (i.e. unoptimized) support+		     for doing alignments - including dynamic adjustment of scores based on sequence quality. +		     Also Blast output parsing.  Partly implemented single linkage clustering, and 		     multiple alignment.  Reading Gene Ontology (GO) annotations (GOA) and 		     definitions\/hierarchy. 		     .@@ -25,7 +26,7 @@  Tested-With:         GHC==6.8.2 Build-Type:          Simple-Build-Depends:       base>3, QuickCheck<2, binary, tagsoup>=0.4, bytestring >= 0.9.1,+Build-Depends:       base>=3 && <4, 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 @@ -40,12 +41,13 @@                      Bio.Sequence.GOA,                      Bio.Sequence.GeneOntology, 		     Bio.Sequence.KEGG,-		     Bio.Sequence.SFF,+		     Bio.Sequence.SFF, Bio.Sequence.SFF_name                      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.Bowtie,                      Bio.Alignment.Soap,                      Bio.Clustering,                      Bio.Util, Bio.Util.Parsex, Bio.Util.TestBase