bio 0.4 → 0.4.4

raw patch · 16 files changed

+595/−101 lines, 16 filesdep ~QuickCheckdep ~basedep ~binaryPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependency ranges changed: QuickCheck, base, binary, tagsoup

API changes (from Hackage documentation)

Files

• Bio/Alignment/AlignData.hs +2/−0
• Bio/Alignment/BED.hs +68/−0
• Bio/Alignment/BlastXML.hs +6/−3
• Bio/Alignment/Matrices.hs +15/−10
• Bio/Alignment/PSL.hs +62/−0
• Bio/Sequence.hs +16/−1
• Bio/Sequence/AminoProperties.hs +133/−0
• Bio/Sequence/Fasta.hs +9/−5
• Bio/Sequence/HashWord.lhs +1/−1
• Bio/Sequence/SFF.hs +140/−38
• Bio/Sequence/SFF_filters.hs +87/−0
• Bio/Sequence/SeqData.hs +27/−13
• Bio/Sequence/TwoBit.hs +2/−3
• Bio/Util/TestBase.hs +2/−11
• README +19/−10
• bio.cabal +6/−6

@@ -41,6 +41,8 @@  -- | Gaps are coded as '*'s, this function removes them, and returns --   the sequence along with the list of gap positions.+--   note that gaps are positioned relative to the *gapped* sequence +--   (contrast to stmassembler/Cluster.hs) extractGaps :: SeqData -> (SeqData,Gaps) extractGaps str = (BC.filter (/='*') str,BC.elemIndices '*' str) 

@@ -0,0 +1,68 @@+-- | Model the BED format, according to the spec at+-- http://genome.ucsc.edu/FAQ/FAQformat#format1++module Bio.Alignment.BED ( BED(..),Dir(..)+                         , readBED, writeBED+                         ) where++import Bio.Sequence.SeqData (Offset)+import Data.Word+import qualified Data.ByteString.Lazy.Char8 as B+import Data.ByteString.Lazy.Char8 (ByteString)++-- | The BED data type Note that the specification allows a variable number of fields, with+-- only the three first required.  This definition requires all fields to be present.+data BED = BED { chrom                :: ByteString+               , chromStart, chromEnd :: Offset+               , name                 :: ByteString+               , score                :: Int                 -- ^ Range 0..1000+               , strand               :: Dir+               , thickStart, thickEnd :: Offset+               , itemRGB              :: (Word8,Word8,Word8) -- ^ Available BED files appear to not+                                                             -- support this format.  RGB is therefore+                                                             -- ignored (read and written as '0')+               , blockSizeStart :: [(Offset,Offset)]         -- ^ Lists of lenght blockCount, blockStarts+                                                             --  are relative to chromStart+               }++-- | Yet another direction data structure.+data Dir = Fwd | Rev deriving Eq++instance Show Dir where+    show Fwd = "+"+    show Rev = "-"++instance Read Dir where+    readsPrec _ ('+':rest) = [(Fwd,rest)]+    readsPrec _ ('-':rest) = [(Rev,rest)]+    readsPrec _ x = error ("Can't parse '"++x++"' as a Dir")++readBED :: FilePath -> IO [BED]+readBED f = (map (unpack1 . myWords) . B.lines) `fmap` B.readFile f+    where unpack1 [c,cs,ce,nm,sc,str,ts,te,_rgb,_bc,bsz,bst] =+              BED c (i cs) (i ce) nm (i sc) (i str) (i ts) (i te) (0,0,0) (zip (i' bsz) (i' bst))+          unpack1 _x = error ("incorrect number of fields in BED record ("++show (length _x)++"):\n"++show _x)++          i :: Read a => ByteString -> a+          i = read . B.unpack+          i' :: Read a => ByteString -> [a]+          i' = read . (++"]") . ('[':) . B.unpack++          myWords = map rmtabs . B.groupBy (const (/='\t'))+          rmtabs x = if not (B.null x) && B.head x == '\t' then B.drop 1 x else x++instance Show BED where+    show = B.unpack . pack1++writeBED :: FilePath -> [BED] -> IO ()+writeBED f = B.writeFile f . B.unlines . map pack1++pack1 :: BED -> ByteString+pack1 (BED c cs ce nm sc str ts te _rgb bszst) =+          let bl = length bszst+              (bsz,bst) = unzip bszst+              b :: Show a => a -> ByteString+              b = B.pack . show+              b':: Show a => [a] -> ByteString+              b'= B.pack . drop 1 . init . show+          in B.intercalate (B.pack "\t") [c, b cs, b ce, nm, b sc, b str, b ts, b te, b (0::Int), b bl, b' bsz, b' bst]

@@ -91,7 +91,10 @@                , h_from = read $ get "Hsp_hit-from"                , h_to   = read $ get "Hsp_hit-to"                , identity = (read $ get "Hsp_identity", read $ get "Hsp_align-len")-               -- blastx has query-frame, tblastn has hit-frame, blastn has both hit and query+               -- blastx has query-frame ±1..3 +               -- tblastn has hit-frame+               -- blastn has both hit and query+               -- tblastx has query-frame = 1, hit-frame ±1..3                , aux = case sections (isTagOpenName "Hsp_hit-frame") ms of                          [] -> mkFrame $ get "Hsp_query-frame"                          [(_o:TagText hf:_c)] -> case sections (isTagOpenName "Hsp_query-frame") ms of @@ -103,8 +106,8 @@     where get = getFrom ms           mkFrame f = Frame (strand' $ signum $ read f) (abs $ read f)           mkStrands h q = Strands (strand' $ read h) (strand' $ read q)+          -- ignore frame also for tblastx hits (it can be reconstructed from location)           strand' :: Int -> Strand-          strand' s = case s of 1 -> Plus; -1 -> Minus-                                _ -> error ("Strand must be +1 or -1, but was: "++show s)+          strand' s = if s > 0 then Plus else Minus  

@@ -34,10 +34,15 @@  import Bio.Alignment.AlignData (Chr) import qualified Data.Map as M+import Data.Char (ord) +genMatrix :: [((Char,Char),Int)] -> M.Map (Chr,Chr) Int+genMatrix = M.fromList . map toW8+    where toW8 ((x,y),i) = ((fromIntegral $ ord x, fromIntegral $ ord y),i)+ -- | BLOSUM45 matrix, suitable for distantly related sequences-blosum45 :: (Char,Char) -> Int-blosum45 m = M.findWithDefault (-5) m $ M.fromList [(('A','A'),5),+blosum45 :: (Chr,Chr) -> Int+blosum45 m = M.findWithDefault (-5) m $ genMatrix [(('A','A'),5),     (('A','B'),-1), (('A','C'),-1), (('A','D'),-2), (('A','E'),-1),     (('A','F'),-2), (('A','G'),0), (('A','H'),-2), (('A','I'),-1),     (('A','K'),-1), (('A','L'),-1), (('A','M'),-1), (('A','N'),-1),@@ -172,8 +177,8 @@     (('Z','W'),-2), (('Z','X'),-1), (('Z','Y'),-2), (('Z','Z'),4)]  -- | The standard BLOSUM62 matrix.-blosum62 :: (Char,Char) -> Int-blosum62 m = M.findWithDefault (-4) m $ M.fromList +blosum62 :: (Chr,Chr) -> Int+blosum62 m = M.findWithDefault (-4) m $ genMatrix      [(('A','A'),4), (('A','B'),-2), (('A','C'),0), (('A','D'),-2),      (('A','E'),-1), (('A','F'),-2), (('A','G'),0), (('A','H'),-2),      (('A','I'),-1), (('A','K'),-1), (('A','L'),-1), (('A','M'),-1),@@ -309,8 +314,8 @@      (('Z','Z'),4)]  -- | BLOSUM80 matrix, suitable for closely related sequences. -blosum80 :: (Char,Char) -> Int-blosum80 m = M.findWithDefault (-6) m $ M.fromList [(('A','A'),5),+blosum80 :: (Chr,Chr) -> Int+blosum80 m = M.findWithDefault (-6) m $ genMatrix [(('A','A'),5),     (('A','B'),-2), (('A','C'),-1), (('A','D'),-2), (('A','E'),-1),     (('A','F'),-3), (('A','G'),0), (('A','H'),-2), (('A','I'),-2),     (('A','K'),-1), (('A','L'),-2), (('A','M'),-1), (('A','N'),-2),@@ -445,8 +450,8 @@     (('Z','W'),-4), (('Z','X'),-1), (('Z','Y'),-3), (('Z','Z'),4)]  -- | The standard PAM30 matrix-pam30 :: (Char,Char) -> Int-pam30 m = M.findWithDefault (-17) m $ M.fromList [(('A','A'),6),+pam30 :: (Chr,Chr) -> Int+pam30 m = M.findWithDefault (-17) m $ genMatrix [(('A','A'),6),     (('A','B'),-3), (('A','C'),-6), (('A','D'),-3), (('A','E'),-2),     (('A','F'),-8), (('A','G'),-2), (('A','H'),-7), (('A','I'),-5),     (('A','K'),-7), (('A','L'),-6), (('A','M'),-5), (('A','N'),-4),@@ -582,8 +587,8 @@     (('Z','Z'),6)]  -- | The standard PAM70 matrix.-pam70 :: (Char,Char) -> Int-pam70 m = M.findWithDefault (-11) m $ M.fromList [(('A','A'),5),+pam70 :: (Chr,Chr) -> Int+pam70 m = M.findWithDefault (-11) m $ genMatrix [(('A','A'),5),     (('A','B'),-1), (('A','C'),-4), (('A','D'),-1), (('A','E'),-1),     (('A','F'),-6), (('A','G'),0), (('A','H'),-4), (('A','I'),-2),     (('A','K'),-4), (('A','L'),-4), (('A','M'),-3), (('A','N'),-2),

@@ -0,0 +1,62 @@+{-| This models the PSL format used by e.g. the alignment tool BLAT.  +    It is a simple, textual representation of (spliced) alignments,+    with tab-separated fields.++    See http://genome.ucsc.edu/FAQ/FAQformat#format2 for details.+-}+module Bio.Alignment.PSL where++import qualified Data.ByteString.Lazy.Char8 as B+import Data.ByteString.Lazy.Char8 (ByteString)++-- encodes the fields in a Blat PSL record+-- NB! laziness issues, if not fully evaluated, this will hang onto the input.+data PSL = PSL   { match, mismatch, repmatch, ncount+                 , qgapcount, qgaplength, tgapcount, tgaplength :: Int +                 , strand :: ByteString+                 , qname :: ByteString, qsize, qstart, qend :: Int+                 , tname :: ByteString, tsize, tstart, tend :: Int+                 , blockcount :: Int, blocksizes, qstarts, tstarts :: [Int]+                 } deriving Show+{-+sread :: ByteString -> Strand+sread b = case B.uncons b of +  Just ('-',_) -> Minus+  Just ('+',_) -> Plus+  _ -> error ("Cant't parse '"++B.unpack b++"' as strand")++swrite :: Strand -> ByteString+swrite Minus = B.pack "-"+swrite Plus  = B.pack "+"+-}++readPSL :: FilePath -> IO [PSL]+readPSL f = parsePSL `fmap` B.readFile f++writePSL :: FilePath -> [PSL] -> IO ()+writePSL f =  B.writeFile f . unparsePSL++parsePSL :: ByteString -> [PSL]+parsePSL s = map parseLine $ B.lines $ dropHeader+  where+    -- dropHeader requires strict adherence to BLAT output.  Perhaps a looser check is better?+    dropHeader | pslHeader `B.isPrefixOf` s = B.drop (B.length pslHeader) s+               | otherwise = error "PSL header mismatch: (Todo: insert specifics here)"+    parseLine l = let fs = B.split '\t' l+                      ri :: Int -> Int+                      ri i = maybe (error ("Can't parse field "++show i++" '"++B.unpack (fs!!i)++"' as an integer")) fst $ B.readInt (fs!!i)+                      rl j = map (\w -> maybe undefined fst (B.readInt w)) $ init $ B.split ',' (fs!!j)+                  in PSL (ri 0) (ri 1) (ri 2) (ri 3) (ri 4) (ri 5) (ri 6) (ri 7) (fs!!8) (fs!!9) (ri 10) (ri 11) (ri 12) (fs!!13) (ri 14) (ri 15) (ri 16) (ri 17) (rl 18) (rl 19) (rl 20)+                           +++unparsePSL :: [PSL] -> ByteString+unparsePSL bs = addHeader $ B.unlines $ map format1 bs+  where addHeader = B.append pslHeader+        format1 :: PSL -> ByteString+        format1 = undefined++pslHeader = B.pack ("psLayout version 3\n\n"+                    ++"match\tmis- \trep. \tN's\tQ gap\tQ gap\tT gap\tT gap\tstrand\tQ        \tQ   \tQ    \tQ  \tT        \tT   \tT    \tT  \tblock\tblockSizes \tqStarts\t tStarts\n"+                    ++"     \tmatch\tmatch\t   \tcount\tbases\tcount\tbases\t      \tname     \tsize\tstart\tend\tname     \tsize\tstart\tend\tcount\n"+                    ++"---------------------------------------------------------------------------------------------------------------------------------------------------------------\n")

@@ -18,8 +18,13 @@     , compl, revcompl, revcompl', Nuc, castToNuc     -- ** Protein sequence functionality     , Amino(..), translate, fromIUPAC, toIUPAC, castToAmino+    -- ** Other utility functions+    , defragSeq, seqmap -    -- * File formats+    -- * File IO+    -- ** Generic sequence reading+    , readNuc, readProt+     -- ** The Fasta file format ("Bio.Sequence.Fasta")     , readFasta, hReadFasta      , writeFasta, hWriteFasta@@ -65,3 +70,13 @@ -- sequence-oriented stuff import Bio.Sequence.Entropy import Bio.Sequence.HashWord++-- | Read nucleotide sequences in any format - Fasta, SFF, FastQ, 2bit, PHD...+readNuc :: FilePath -> IO [Sequence Nuc]+readNuc = undefined+          -- check file contents+          -- magic number++-- | Read protein sequences in any supported format (i.e. Fasta)+readProt :: FilePath -> IO [Sequence Amino]+readProt xs = map castSeq `fmap` readFasta xs

@@ -0,0 +1,133 @@+{-| Encodes a table of amino acid properties. +    Based on Livingstone & Barton, CABIOS, 9, 745-756, 1993, as seen at:+    http://www.compbio.dundee.ac.uk/user/ws-dev1/jalview/latest/help/html/misc/aaproperties.html+   NB: based on the graphic, not the table (in which P is polar, but T is not)+|-}++module Bio.Sequence.AminoProperties where++import Prelude hiding (or)+import Bio.Sequence.SeqData++type AAProp = Amino -> Bool++oneOf = flip elem+or f g = \x -> f x || g x+  +aliphatic, aromatic, hydrophobic, polar, small, tiny, charged, negative, positive :: AAProp++-- NB: Cysteine has two variants: C s-s (small, non-polar) and C s-h (tiny, polar)+-- http://www.russell.embl-heidelberg.de/aas/aas.html++aliphatic   = oneOf [Ile,Leu,Val]+aromatic    = oneOf [Phe,Tyr,Trp,His]+hydrophobic = aromatic `or` aliphatic `or` oneOf [Cys,Ala,Gly,Thr,Met,Lys]+polar       = charged `or` oneOf [Asn,Ser,Thr,Tyr,Trp,Gln] -- and C s-h+small       = tiny `or` oneOf [Pro,Cys,Val,Thr,Asp,Asn] -- +Asx?+tiny        = oneOf [Ala,Gly,Ser] -- C s-h+charged     = negative `or` positive+negative    = oneOf [Asp,Glu]+positive    = oneOf [Lys,His,Arg]++-- Based on Kyte and Doolittle, according to Wikipedia+hydropathy :: Amino -> Double+hydropathy x = case x of+    Ala ->  1.8 +    Arg -> -4.5 +    Asn -> -3.5 +    Asp -> -3.5 +    Cys ->  2.5 +    Gln -> -3.5 +    Glu -> -3.5 +    Gly -> -0.4 +    His -> -3.2 +    Ile ->  4.5 +    Leu ->  3.8 +    Lys -> -3.9 +    Met ->  1.9 +    Phe ->  2.8 +    Pro -> -1.6 +    Ser -> -0.8 +    Thr -> -0.7 +    Trp -> -0.9 +    Tyr -> -1.3 +    Val ->  4.2+    _   ->  0+    +mass :: Amino -> Double+mass x = case x of +  -- tiny+  Gly ->  57.021464+  Ala ->  71.037114 +  Ser ->  87.032029 +  -- small+  Pro ->  97.052764 +  Val ->  99.068414+  Thr -> 101.04768 +  Cys -> 103.00919 +  Asn -> 114.04293 +  Asp -> 115.02694+  +  Leu -> 113.08406 +  Ile -> 113.08406 +  Gln -> 128.05858+  Lys -> 128.09496+  Glu -> 129.04259+  Met -> 131.04048+  His -> 137.05891+  Phe -> 147.06841+  Arg -> 156.10111+  Tyr -> 163.06333+  Trp -> 186.07931+  _   -> 0++-- | The propensities for forming secondary structures+--  From Zvelebil and Baum: Understanding Bioinformatics, Chapter 11+--  citing Chou and Fasman.+--  Today, more complex methods like GOR are recommended instead.+helixP, strandP :: Amino -> Double+helixP x = case x of+  Ala -> 1.42+  Cys -> 0.70+  Asp -> 1.01+  Glu -> 1.51+  Phe -> 1.13+  Gly -> 0.61+  His -> 1.00+  Ile -> 1.08+  Lys -> 1.16+  Leu -> 1.21+  Met -> 1.45+  Asn -> 0.67+  Pro -> 0.57+  Gln -> 1.11+  Arg -> 0.98+  Ser -> 0.77+  Thr -> 0.83+  Val -> 1.06+  Trp -> 1.08+  Tyr -> 0.69+  _   -> 0    -- sensible?+  +strandP x = case x of  +  Ala -> 0.83+  Cys -> 1.19+  Asp -> 0.54+  Glu -> 0.37+  Phe -> 1.38+  Gly -> 0.75+  His -> 0.87+  Ile -> 1.60+  Lys -> 0.74+  Leu -> 1.30+  Met -> 1.05+  Asn -> 0.89+  Pro -> 0.55+  Gln -> 1.10+  Arg -> 0.93+  Ser -> 0.75+  Thr -> 1.19+  Val -> 1.70+  Trp -> 1.37+  Tyr -> 1.40+  _   -> 0

@@ -162,11 +162,15 @@  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)))-          f [] = error "mkQual: empty quality data"-          int = fromIntegral . fst . fromJust' . B.readInt-          fromJust' = maybe (error "Error in qual format") id+    where f [] = error "mkQual on empty input - this is impossible"+          f (l:ls) = Seq (B.drop 1 l) B.empty +                     (Just $ BB.pack $ go 0 ls)+          isDigit x = x <= 58 && x >= 48 +          go i (s:ss) = case BB.uncons s of Just (c,rs) -> if isDigit c then go (c - 48 + 10*i) (rs:ss)+                                                           else let rs' = BB.dropWhile (not . isDigit) rs+                                                                in if BB.null rs' then i : go 0 ss else i : go 0 (rs':ss)+                                            Nothing -> i : go 0 ss+          go _ [] = []  -- | Split lines into blocks starting with '>' characters --   Filter out # comments (but not semicolons?)

@@ -192,7 +192,7 @@ k2n :: Integral k => Int -> k -> SeqData k2n k = fromStr . k2n' k -k2n' k i = if k==1 then [unval i] +k2n' k i = if k==1 then [unval i]            else let (q,r) = i `divMod` (4^(k-1)) in unval q : k2n' (k-1) r  {-

@@ -13,9 +13,11 @@  module Bio.Sequence.SFF ( SFF(..), CommonHeader(..)                         , ReadHeader(..), ReadBlock(..)-                        , readSFF, writeSFF, writeSFF'-                        , sffToSequence-                        , test, convert+                        , readSFF, writeSFF, writeSFF', recoverSFF+                        , sffToSequence, trim, trimFromTo, trimKey+                        , baseToFlowPos, flowToBasePos+                        , test, convert, flowgram+                        , packFlows, unpackFlows                         , Flow, Qual, Index, SeqData, QualData                         , ReadName (..), decodeReadName, encodeReadName                         ) where@@ -29,14 +31,14 @@ import qualified Data.ByteString as B import qualified Data.ByteString.Char8 as BC import Data.ByteString (ByteString)-import Control.Monad (when,replicateM,replicateM_)+import Control.Monad (when,replicateM,replicateM_,liftM)  import Data.List (intersperse) import Data.Binary import Data.Binary.Get (getByteString,getLazyByteString) import qualified Data.Binary.Get as G import Data.Binary.Put (putByteString,putLazyByteString)-+import Data.Char (toUpper, toLower) import Text.Printf (printf) import System.IO @@ -51,42 +53,79 @@ versions :: [Int32] versions = [1] - readSFF :: FilePath -> IO SFF readSFF f = return . decode =<< LB.readFile f +trimKey :: CommonHeader -> Sequence Nuc -> Maybe (Sequence Nuc)+trimKey ch (Seq n s q) = let (k,s2) = LB.splitAt (fromIntegral $ key_length ch) s+                          in if LBC.map toLower k==LBC.map toLower (LB.fromChunks [key ch]) +                             then Just $ Seq n s2 (liftM (LB.drop (fromIntegral $ key_length ch)) q)+                             else Nothing -- error ("Couldn't match key in sequence "++LBC.unpack n++" ("++LBC.unpack k++" vs. "++BC.unpack (key ch)++")!")+ 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) = 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)++")!")+    where r2s r = clip (read_header r, bases r, quality r)+          clip (h, s, q) = let (left,right) = (clip_qual_left h,clip_qual_right h)+                               split x = let (a,b) = LB.splitAt (fromIntegral right) x +                                             (c,d) = LB.splitAt (fromIntegral left-1) a+                                         in [c,d,b]+                           in {- trim_key $ -} Seq (LB.fromChunks [read_name h ,BC.pack (" qclip: "++show left++".."++show right)])+                                  (let [a,b,c] = split s in LBC.concat [LBC.map toLower a,LBC.map toUpper b,LBC.map toLower c])+                                  (Just q) +-- trimming the flowgram is necessary, but how to deal with the shift in flow+-- sequence - i.e. what to do when trimming "splits" a flow into trimmed/untrimmed bases?+-- | Trim a read to specific sequence position.+-- The current implementation has the unintended side effect of always trimming the flowgram down to a basecalled position.+trimFromTo :: (Integral i) => i -> i -> ReadBlock -> ReadBlock+trimFromTo l r rd = let trim_seq = LB.drop (fromIntegral l) . LB.take (fromIntegral r)+                        trim_seq' = B.drop (fromIntegral l) . B.take (fromIntegral r)+                        trim_flw = B.drop ((2*) $ fromIntegral $ baseToFlowPos rd l) . B.take ((2*) $ fromIntegral $ baseToFlowPos rd r)+                    in rd { read_header = read_header rd -- FIXME: Update num_bases?+                          , flow_data = trim_flw (flow_data rd)+                          , flow_index = trim_seq' (flow_index rd)+                          , bases = trim_seq (bases rd)+                          , quality = trim_seq (quality rd)+                          }++-- | Trim a read according to clipping information+trim :: ReadBlock -> ReadBlock+trim rb = let rh = read_header rb in trimFromTo (clip_qual_left rh-1) (clip_qual_right rh) rb++-- | Convert a flow position to the corresponding sequence position+flowToBasePos :: Integral i => ReadBlock -> i -> Int+flowToBasePos rd fp = length $ takeWhile (<fp) $ scanl (+) 0 $ map fromIntegral $ B.unpack $ flow_index rd++-- | Convert a sequence position to the corresponding flow position+baseToFlowPos :: Integral i => ReadBlock -> i -> Int+baseToFlowPos rd sp = sum $ map fromIntegral $ B.unpack $ B.take (fromIntegral sp) $ flow_index rd++recoverSFF :: FilePath -> IO SFF+recoverSFF f = return . unRecovered . decode =<< LB.readFile f+ -- | 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 ()+--   stream of 'ReadBlock's.  Returns the number of reads written.+writeSFF' :: FilePath -> SFF -> IO Int 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+  return $ fromIntegral c  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)+  return $! (c+1)  data RBI = RBI Int ReadBlock @@ -131,20 +170,11 @@     get = do       -- Parse CommonHeader       chead <- get-      -- Parse ReadHeaders+      -- Get the ReadBlocks       rds <- replicateM (fromIntegral (num_reads chead))                                    (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 <- 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)+                                      getRB chead rh                                    )       return (SFF chead rds) @@ -152,14 +182,28 @@       put hd       mapM_ (put . RBI (fromIntegral $ flow_length hd)) rds +{-# INLINE getRB #-}+getRB :: CommonHeader -> ReadHeader -> Get ReadBlock+getRB chead rh = do+  let nb = fromIntegral $ num_bases rh+      nb' = fromIntegral $ num_bases rh+      fl = fromIntegral $ flow_length chead+  fg <- getByteString (2*fl)+  fi <- 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 fg fi bs qty)+ -- | A ReadBlock can't be an instance of Binary directly, since it depends on --   information from the CommonHeader. putRB :: Int -> ReadBlock -> Put putRB fl rb = do   put (read_header rb)-  mapM_ put (flowgram rb)+  putByteString (flow_data rb)   -- ensure that flowgram has correct lenght-  replicateM (2*(fl-length (flowgram rb))) (put (0::Word8))+  replicateM (2*fl-B.length (flow_data rb)) (put (0::Word8))   putByteString (flow_index rb)   putLazyByteString (bases rb)   putLazyByteString (quality rb)@@ -167,13 +211,19 @@       l = (fl*2+nb*3) `mod` 8   when (l > 0) (pad (8-l)) --- | What the name and type says.-decodeArray :: ByteString -> [Flow]-decodeArray = dec . map fromIntegral . B.unpack +-- | Unpack the flow_data field into a list of flow values+unpackFlows :: ByteString -> [Flow]+unpackFlows = dec . map fromIntegral . B.unpack      where dec (d1:d2:rest) = d1*256+d2 : dec rest           dec [] = []           dec _  = error "odd flowgram length?!" +-- | Pack a list of flows into the corresponding binary structure (the flow_data field)+packFlows :: [Flow] -> ByteString+packFlows = B.pack . map fromIntegral . merge +  where merge (x:xs) = let (a,b) = x `divMod` 256 in a:b:merge xs+        merge [] = []+ -- ---------------------------------------------------------- -- | SFF has a 31-byte common header --   Todo: remove items that are derivable (counters, magic, etc)@@ -267,20 +317,72 @@ -- ---------------------------------------------------------- -- | This contains the actual flowgram for a single read. data ReadBlock = ReadBlock {-      read_header                :: ReadHeader+      read_header                :: ! ReadHeader     -- The data block-    , flowgram                   :: [Flow]-    , flow_index                 :: ByteString-    , bases                      :: SeqData-    , quality                    :: QualData+    , flow_data                  :: ! ByteString -- nb! use unpackFlows for this+    , flow_index                 :: ! ByteString+    , bases                      :: ! SeqData+    , quality                    :: ! QualData     } +flowgram :: ReadBlock -> [Flow]+flowgram = unpackFlows . flow_data+ instance Show ReadBlock where     show (ReadBlock h f i b q) =         show h ++ unlines (map ("     "++) -            ["flowgram:\t"++show f+            ["flowgram:\t"++show (unpackFlows f)             , "index:\t"++(concat . intersperse " " . map show . B.unpack) i             , "bases:\t"++LBC.unpack b             , "quality:\t"++(concat . intersperse " " . map show . LB.unpack) q             , ""             ])++-- ------------------------------------------------------------+-- | RSFF wraps an SFF to provide an instance of Binary with some more error checking.+data RSFF = RSFF { unRecovered :: SFF }++instance Binary RSFF where +    get = do+      -- Parse CommonHeader+      chead <- get+      -- Get the first read block+      r1 <- do rh <- get +               getRB chead rh+      -- Get subsequent read blocks+      rds <- replicateM (fromIntegral (num_reads chead))+                                   (do rh <- getSaneHeader (take 4 $ BC.unpack $ read_name $ read_header r1)+                                       getRB chead rh)+      return (RSFF $ SFF chead (r1:rds))+    put = error "You should not serialize an RSFF"++-- | This allows us to decode the constant parts of the read header for verifying its correcness.+data PartialReadHeader = PartialReadHeader {+      _pread_header_lenght                    :: Int16 +    , _pname_length                           :: Int16+    , _pnum_bases                             :: Int32+    , _pclip_qual_left, _pclip_qual_right+    , _clip_adapter_left, _pclip_adapter_right :: Int16+    , _pread_name                              :: ByteString -- length four+}++instance Binary PartialReadHeader where+    get = do { rhl <- get; nl <- get; nb <- get; ql <- get; qr <- get; al <- get; ar <- get; rn <- getByteString 4 +             ; return (PartialReadHeader rhl nl nb ql qr al ar rn) }+    put = error "You should not serialize a PartialReadHeader"++-- | Ensure that the header we're decoding matches our expectations.+getSaneHeader :: String -> Get ReadHeader+getSaneHeader prefix = do+  buf <- getLazyByteString 20+  decodeSaneH prefix buf  ++decodeSaneH :: String -> LBC.ByteString -> Get ReadHeader+decodeSaneH prefix buf = do+  let PartialReadHeader rhl nl _nb _ql _qr _al _ar rn = decode buf+  if rhl >= 20 && nl > 0 && all id (zipWith (==) prefix (BC.unpack rn))+      then do buf2 <- getLazyByteString (fromIntegral rhl-20)+              return (decode $ LB.concat [buf,buf2])+      else do x <- getLazyByteString 1 -- error "skip one byte, try again"+              decodeSaneH prefix (LBC.concat [buf,x])+

@@ -0,0 +1,87 @@+-- | This implements a number of filters used in the Titanium pipeline+module Bio.Sequence.SFF_filters where++import Bio.Sequence.SFF (ReadBlock(..), ReadHeader(..), flowToBasePos, flowgram)+import qualified Data.ByteString.Lazy as B+import qualified Data.ByteString.Lazy.Char8 as BL+import Data.List (tails)++-- Ti uses a set of filters, described in the (something) manual.+-- (GS Run Processor Application, section 3.2.2++)++-- ** Discarding filters **++-- | DiscardFilters determine whether a read is to be retained or discarded+type DiscardFilter = ReadBlock -> Bool -- True to retain, False to discard++filter_dots, filter_mixed, filter_key, filter_empty :: DiscardFilter++filter_empty rb = num_bases (read_header rb) >= 5++filter_key rb = ("TCAG"==) $ take 4 $ BL.unpack $ bases rb++-- | 3.2.2.1.2 The "dots" filter discards sequences where the last positive flow is +--   before flow 84, and flows with >5% dots (i.e. three successive noise values) +--   before the last postitive flow. (Interpreted as 5% of called sequence length is Ns?)+filter_dots rb = let dots = BL.length $ BL.filter (=='N') $ bases rb+                 in fromIntegral dots / fromIntegral (BL.length $ bases rb) < (0.05 :: Double)++-- | 3.2.2.1.3 The "mixed" filter discards sequences with more than 70% positive flows.  +--   Also, discard with <30% noise, >20% middle (0.45..0.75) or <30% positive.+filter_mixed rb = let fs = dropWhile (<50) . reverse . flowgram $ rb+                      fl = dlength fs+                  in and [ (dlength (filter (>50) fs) / fl) > 0.7+                         , (dlength (filter (<45) fs) / fl) > 0.3 +                         , (dlength (filter (>75) fs) / fl) > 0.3+                         , (dlength (filter (\f -> f<=75 && f>=45) fs) / fl) < 0.2+                         ]+-- | Discard a read if the number of untrimmed flows is less than n (n=186 for Titanium)+filter_length :: Int -> DiscardFilter+filter_length n rb = length (flowgram rb) >= n++-- ** Trimming filters **++-- | TrimFilters modify the read, typically trimming it for quality+type TrimFilter = ReadBlock -> ReadBlock+filter_sigint, filter_qual20 :: TrimFilter++-- | 3.2.2.1.4 Signal intensity trim - trim back until <3% borderline flows (0.5..0.7).+--   Then trim borderline values or dots from the end (use a window).+filter_sigint rb = clipFlows rb (sigint rb)++sigint :: ReadBlock -> Int+sigint rb = let bs = scanl (\(n,m,_) f -> if f >= 50 && f <= 70 then (n+1,m+1,f) else (n,m+1,f)) (0,0,0) $ flowgram rb +            in length $ reverse +                   $ dropWhile (\(_,_,f) -> f<=70)+                   $ dropWhile (\(n,m,_)->n `div` (m*1000) > (30::Int)) $ reverse bs+++-- 3.2.2.1.5 Primer filter and 3.2.2.1.6 Trimback valley filter are ignored, +-- since we don't know how to detect the primer, and don't understand the description of tbv.++-- | 3.2.2.1.7 Quality score trimming trims using a 10-base window until a Q20 average is found. +filter_qual20 rs = clipSeq rs $ qual20 rs++qual20 :: ReadBlock -> Int+qual20 rs = (fromIntegral $ num_bases $ read_header rs) +            - (length . takeWhile (<20) . map (avg . take 10) . tails . reverse . B.unpack $ quality rs)++-- ** Utility functions **++-- | List length as a double (eliminates many instances of fromIntegral)+dlength :: [a] -> Double+dlength = fromIntegral . length++-- | Calculate average of a list+avg :: Integral a => [a] -> Double+avg xs = fromIntegral (sum xs) / dlength xs++-- | Translate a number of flows to position in sequence, and update clipping data accordingly+clipFlows :: ReadBlock -> Int -> ReadBlock+clipFlows rb n = clipSeq rb (flowToBasePos rb n)++-- | Update clip_qual_right if more severe than previous value+clipSeq :: ReadBlock -> Int -> ReadBlock+clipSeq rb n' = let n = fromIntegral n' +                    rh = read_header rb +                in if clip_qual_right rh <= n then rb else rb { read_header = rh {clip_qual_right = n }}

@@ -30,7 +30,7 @@     , fromStr, toStr        -- * Sequence utilities-    , defragSeq+    , defragSeq, seqmap, castSeq        -- * Nucleotide functionality       -- | Nucleotide sequences contain the alphabet [A,C,G,T].@@ -55,9 +55,9 @@ 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) @@ -83,13 +83,15 @@ data Nuc data Unknown --- | Phantom type functionality+-- | Phantom type functionality, unchecked conversion between sequence types+castSeq :: Sequence a -> Sequence b+castSeq (Seq h d mq) = Seq h d mq+ castToNuc :: Sequence a -> Sequence Nuc-castToNuc (Seq h d mq) = Seq h d mq+castToNuc = castSeq  castToAmino :: Sequence a -> Sequence Amino-castToAmino (Seq h d mq) = Seq h d mq-+castToAmino = castSeq  -- | A more arranged show instance for Sequences reassembling the display of --   the fasta-format@@ -228,7 +230,7 @@ dropSkip :: (Eq a) => [a] -> Int -> [[a]] -> [a] dropSkip []         _ _     = [] dropSkip xs         0 _     = xs-dropSkip str@(x:xs) n skips =+dropSkip str@(_:xs) n skips =     let n' = case (testPrefixes str skips) of                Nothing  -> n-1                Just len -> n + len - 1@@ -239,12 +241,12 @@ -- 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 = +insert xs         0   e    _     = e ++ xs+insert str@(x:xs) pos e skips =      let pos' = case (testPrefixes str skips) of                Nothing  -> pos-1                Just len -> pos + len - 1-    in x : insert xs pos' elem skips+    in x : insert xs pos' e skips          -- Test, if one of the prefixes is the prefix of xs. If yes, return it's length -- else Nothing  @@ -319,9 +321,21 @@ -- 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+    where defragB = B.fromChunks . (: []) . BBS.concat . B.toChunks           defragBB = BB.fromChunks . (: []) . BBS.concat . BB.toChunks ++-- | map over sequences, treating them as a sequence of (char,word8) pairs.+--   This will work on sequences without quality, as long as the function doesn't +--   try to examine it.+--   The current implementation is not very efficient.+seqmap :: ((Char,Qual) -> (Char,Qual)) -> Sequence t -> Sequence t+seqmap f (Seq n ds Nothing) = Seq n xs Nothing +    where xs = B.map (\x -> fst $ f (x,e)) ds+          e = error "seqmap: tried to access qual data, but none available"+seqmap f (Seq n ds (Just qs)) = Seq n (B.pack xs) (Just $ BB.pack ys)+    where (xs,ys) = unzip $ map f $ zip (B.unpack ds) (BB.unpack qs)+ ------------------------------------------------------------ -- Nucleotide (DNA, RNA) specific stuff ------------------------------------------------------------@@ -335,7 +349,7 @@  -- | Calculate the reverse complent for SeqData only. revcompl' :: SeqData -> SeqData-revcompl' = B.reverse . B.map compl+revcompl' = B.map compl . B.reverse  -- | Complement a single character.  I.e. identify the nucleotide it  --   can hybridize with.  Note that for multiple nucleotides, you usually@@ -409,7 +423,7 @@  -- | Convert a sequence in IUPAC format to a list of amino acids. fromIUPAC :: SeqData -> [Amino]-fromIUPAC = map (maybe Xaa id . flip lookup iupac') . B.unpack+fromIUPAC = map (maybe Xaa id . flip lookup iupac' . toUpper) . B.unpack  iupac :: [(Amino,Char)] iupac = [(Ala,'A')        ,(Arg,'R')        ,(Asn,'N')

@@ -34,9 +34,8 @@ import Data.List import Data.Bits -import Test.QuickCheck hiding (check)    -- QC 1.0--- import Test.QuickCheck hiding ((.&.)) -- QC 2.0-+-- import Test.QuickCheck hiding (check)    -- QC 1.0+import Test.QuickCheck hiding ((.&.)) -- QC 2.0  -- constants default_magic, default_version :: Word32

@@ -78,16 +78,13 @@  instance Arbitrary Word8 where     arbitrary = choose (0,255)-    coarbitrary _ = id  instance Arbitrary Nucleotide where     arbitrary = elements (map N "aaacccgggtttn")-    coarbitrary _ = id  instance Arbitrary Quality where     arbitrary = do c <- choose (0,60)                    return (Q c)-    coarbitrary _ = id  instance Arbitrary ESTq where     arbitrary = do n <- choose (1,100)@@ -95,18 +92,15 @@                    q <- vector n                    return $ Eq $ Seq (fromStr "qctest")                               (fromStr $ map fromN s) (Just $ pack $ map fromQ q)-    coarbitrary _ = id  instance Arbitrary EST where     arbitrary = do n <- choose (1,100)                    s <- vector n                    return $ E $ Seq (fromStr "qctest")                               (fromStr $ map fromN s) Nothing-    coarbitrary _ = id -instance Arbitrary Char where-    arbitrary = elements (['A'..'Z']++['a'..'z']++" \t\n\r")-    coarbitrary _ = id+-- instance Arbitrary Char where+--    arbitrary = elements (['A'..'Z']++['a'..'z']++" \t\n\r")  instance Arbitrary EST_short where     arbitrary = do let n = 200@@ -114,7 +108,6 @@                    q <- vector n                    return $ ES $ Seq (fromStr "qctest")                               (fromStr $ map fromN s) (Just $ pack $ map fromQ q)-    coarbitrary _ = id  instance Arbitrary EST_long where     arbitrary = do let n = 1000@@ -122,11 +115,9 @@                    q <- vector n                    return $ EL $ Seq (fromStr "qctest")                               (fromStr $ map fromN s) (Just $ pack $ map fromQ q)-    coarbitrary _ = id  -- 1000 ESTs of length 1000 instance Arbitrary EST_set where     arbitrary = do let n = 1000                    s <- vector n                    return (ESet $ map (\(EL x) -> x) s)-    coarbitrary _ = id

@@ -7,26 +7,35 @@ Current list of features includes: a Sequence data type supporting protein and nucleotide sequences and conversion between them, 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-multiple alignment.+phd formats, and also FastQ.  Rudimentary support for doing alignments+- including dynamic adjustment of scores based on sequence quality -+and Blast output parsing.  Partly implemented single linkage+clustering, and multiple alignment.  Reading and writing 454-style SFF+sequences, with a bunch of useful (and not-so-useful) operations on+them.  To install, you need to acquire a working GHC (possibly other Haskell system).  You also need the following external libraries: -  QuickCheck   - for unit tests+  QuickCheck   - for unit tests ('make test' to run them)   binary       - mainly for dealing with the TwoBit sequence format   tagsoup      - for parsing XML output from Blast+  parsec       - for parsing ACE, Bowtie and Soap output  You should be able to get what you need from <http://hackage.haskell.org/>. -You can then build with 'make', doing either 'make install' if you can sudo, or 'make user_install' if you can not.  Of course, the Makefile just proxies for-the regular Cabal routine, which will work just as well:+The easiest way these days, is to use 'cabal'.  You can get away with -    runhaskell Setup configure-    runhaskell Setup build-    sudo runhaskell Setup install+  cabal install bio++As an alternative, you can build with 'make', doing either 'make+install' if you can sudo, or 'make user_install' if you can not.  Of+course, the Makefile just proxies for the regular Cabal routine, which+will work just as well:++  runhaskell Setup configure+  runhaskell Setup build+  sudo runhaskell Setup install  (Use --prefix=$HOME and remove the sudo, if you don't want to install as root.) 

@@ -1,5 +1,5 @@ Name:                bio-Version:             0.4+Version:             0.4.4 License:             LGPL License-file:        LICENSE Author:              Ketil Malde@@ -24,13 +24,11 @@                      The Darcs repository is at: <http://malde.org/~ketil/biohaskell/biolib>. Homepage:            http://blog.malde.org/index.php/the-haskell-bioinformatics-library/ -Tested-With:         GHC==6.8.2+Tested-With:         GHC==6.12.1 Build-Type:          Simple-Build-Depends:       base>=3 && <4, QuickCheck<2, binary, tagsoup>=0.4, bytestring >= 0.9.1,+Build-Depends:       base>=4 && <5, QuickCheck>=2, binary >=0.4 && <0.5, tagsoup>=0.4 && <0.8, bytestring >= 0.9.1,                      containers, array, parallel, parsec, random, old-time, mtl--- add fps for ghc 6.4.2; change imports in Bio/Sequence/TwoBit.hs if you want QC 2 --- We omit the debian/ and Test/ files because those are for development, not installation. Data-Files:          README  Exposed-modules:     Bio.Sequence,@@ -41,7 +39,8 @@                      Bio.Sequence.GOA,                      Bio.Sequence.GeneOntology, 		     Bio.Sequence.KEGG,-		     Bio.Sequence.SFF, Bio.Sequence.SFF_name+		     Bio.Sequence.AminoProperties,+		     Bio.Sequence.SFF, Bio.Sequence.SFF_name, Bio.Sequence.SFF_filters                      Bio.Alignment.BlastData, Bio.Alignment.BlastFlat,                      Bio.Alignment.Blast, Bio.Alignment.BlastXML,                      Bio.Alignment.AlignData, Bio.Alignment.Matrices,@@ -49,6 +48,7 @@                      Bio.Alignment.Multiple, Bio.Alignment.ACE,                      Bio.Alignment.Bowtie,                      Bio.Alignment.Soap,+		     Bio.Alignment.BED, Bio.Alignment.PSL                      Bio.Clustering,                      Bio.Util, Bio.Util.Parsex, Bio.Util.TestBase 		 Bio.Location.Strand, Bio.Location.Position,