bio 0.4.8 → 0.5

raw patch · 20 files changed

+1520/−103 lines, 20 filesdep +processdep ~basedep ~bytestringnew-component:exe:fastoutnew-component:exe:flxnew-component:exe:frecovernew-component:exe:frenamenew-component:exe:qcPVP ok

version bump matches the API change (PVP)

Dependencies added: process

Dependency ranges changed: base, bytestring

API changes (from Hackage documentation)

Files

• Bio/Alignment/Test.hs +362/−0
• Bio/Clustering/Test.hs +38/−0
• Bio/GFF3/Test.hs +64/−0
• Bio/Location/Test.hs +326/−0
• Bio/Sequence.hs +3/−2
• Bio/Sequence/Fasta.hs +6/−11
• Bio/Sequence/SFF.hs +69/−30
• Bio/Sequence/SFF_filters.hs +115/−35
• Bio/Sequence/SeqData.hs +4/−0
• Bio/Sequence/Test.hs +119/−0
• Bio/Sequence/TwoBit.hs +1/−1
• Bio/Util/Test.hs +19/−0
• Bio/Util/TestBase.hs +117/−0
• Makefile +61/−0
• Test.hs +38/−0
• bio.cabal +96/−24
• examples/FRecover.hs +13/−0
• examples/FRename.hs +29/−0
• examples/FastOut.hs +19/−0
• examples/Flx.hs +21/−0

@@ -0,0 +1,362 @@+-- Tests for alignments++module Bio.Alignment.Test where+import Control.Exception (bracket)+import Control.Monad+import Data.Char (isSpace)+import Data.Ord (comparing)+import System.Cmd+import System.Directory+import System.IO+import System.IO.Unsafe++import qualified Bio.Location.ContigLocation as CLoc+import Bio.Location.OnSeq+import qualified Bio.Location.SeqLocation as SeqLoc+import Bio.Location.Strand+import Bio.Sequence+import Bio.Util.TestBase+import Bio.Alignment.AlignData (showalign)+import Bio.Alignment.Matrices as M++import Test.QuickCheck+import Bio.Alignment.QAlign as Q+import Bio.Alignment.AAlign as A++import Bio.Alignment.AlignData (toStrings, extractGaps, insertGaps)+import Bio.Alignment.Multiple++import qualified Bio.Alignment.Soap as Soap++import Data.List (intersperse, sort, sortBy, nub)+import qualified Data.ByteString.Lazy as B+import qualified Data.ByteString.Lazy.Char8 as BC++-- default blastn parameters+mx = Q.qualMx+amx = M.blastn_default+g  = (-5,-2)++tests :: [Test]+--           .........o.........o.........o+tests = [ T "gapped assembly invariant" prop_asm_gaps++        , T "global reverse"        (prop_reverse_g mx g)+        , T "global reverse w/qual" (prop_reverse_g_qual mx g)+        , T "local reverse"         (prop_reverse_l mx g)+        , T "local reverse w/qual"  (prop_reverse_l_qual mx g)+        , T "global recovery"       (prop_recover_g mx g)++        , T "global score <= local" (prop_global_local mx g)+        , T "global score <= overlap" (prop_global_overlap mx g)+        , T "overlap score <= local" (prop_overlap_local mx g)++        , T "overlap score = align" (prop_overlap_score_align mx g)++        , T "s == s => equal scores" (prop_equal_scores mx g)++        -- should give same result without n's(?)+        , T "global aalign = qualign" prop_AAlign_QAlign_g+        , T "local aalign = qualign"  prop_AAlign_QAlign_l++        -- are these really correct?+        , T "decr qual => decr score (l)" (prop_quality_dec_l mx (-100,-100))+        , T "decr qual => decr score (g)" (prop_quality_dec_g mx (-100,-100))++        , T "indirect recover (g)" prop_indirect_recover++        , T "soap parsing" prop_SoapAlign_invertparse+        ]++-- pretty ugly, to ensure the gaps are correctly generated+prop_asm_gaps :: (String,[Int]) -> Bool+prop_asm_gaps (s,gs) = let (str,gaps) = (BC.pack $ filter (/='-') $ filter (/='*') s+                                        , nub $ takeWhile (<= BC.length str) $ dropWhile (<0) $ map fromIntegral $ sort gs)+                       in (str,gaps) == (extractGaps . insertGaps '*' $ (str,gaps))++-- these aren't entirely equivalent, as the different column generators or selector+-- sometimes make different, but equally scoring, choices.  Should be fixed, but in+-- the meantime, checking scores to within an epsilot should be okay.+-- (May fail if two n's are aligned?)+prop_AAlign_QAlign_g (E s1) (E s2) = +    let qmx = qualMx 22 22+        a = (A.global_align qmx g s1 s2)+        q = (Q.global_align mx g s1 s2)+    in if abs (fst a - fst q) < 0.001 then True +       else error ("\n"++show (fst a)++"\n"++ showalign (snd a)+                   ++"\n"++show (fst q) ++ "\n"++ showalign (snd q))++prop_AAlign_QAlign_l (E s1) (E s2) =+    let qmx = qualMx 22 22+        a = (A.local_align qmx g s1 s2)+        q = (Q.local_align mx g s1 s2)+    in if abs (fst a - fst q) < 0.001 then True +       else error ("\n"++show (fst a)++"\n"++ showalign (snd a)+                   ++"\n"++show (fst q) ++ "\n"++ showalign (snd q))++-- Check reverse without quality values+prop_reverse_g mx g (E s1) (E s2) =+    abs (Q.global_score mx g s1 s2 +         - Q.global_score mx g (revcompl s1) (revcompl s2)) < 0.1++prop_reverse_l mx g (E s1) (E s2) =+    abs (Q.local_score mx g s1 s2 +         - Q.local_score mx g (revcompl s1) (revcompl s2)) < 0.1++-- Check reverse with quality+prop_reverse_g_qual mx g (Eq s1) (Eq s2) =+    abs (Q.global_score mx g s1 s2 +           - Q.global_score mx g (revcompl s1) (revcompl s2)) < 0.1++prop_reverse_l_qual mx g (Eq s1) (Eq s2) =+    abs (Q.local_score mx g s1 s2 +           - Q.local_score mx g (revcompl s1) (revcompl s2)) < 0.1++-- Check that alignments are produced from the correct sequences+prop_recover_g mx g (Eq s1) (Eq s2) =+    let (s1',s2') = toStrings $ snd $ Q.global_align mx g s1 s2+    in filter (/='-') s1' == toStr (seqdata s1)+           && filter (/='-') s2' == toStr (seqdata s2)++-- global <= overlap <= local+-- Global score never exceeds optimal local score+prop_global_local mx g (Eq s1) (Eq s2) =+    Q.global_score mx g s1 s2 <= Q.local_score mx g s1 s2++-- Global score never exceeds optimal overlap score+prop_global_overlap mx g (Eq s1) (Eq s2) =+    Q.global_score mx g s1 s2 <= Q.overlap_score mx g s1 s2++-- Overlap score never exceeds optimal local score+prop_overlap_local mx g (Eq s1) (Eq s2) =+    Q.overlap_score mx g s1 s2 <= Q.local_score mx g s1 s2++prop_overlap_score_align mx g (Eq s1) (Eq s2) =+    abs (Q.overlap_score mx g s1 s2 - fst (Q.overlap_align mx g s1 s2)) < 0.1++prop_equal_scores mx g (Eq s1) = +    let qg = Q.global_score mx g s1 s1+        ql = Q.local_score mx g s1 s1+        qo = Q.overlap_score mx g s1 s1+    in qg == ql && ql == qo++-- Sinking score with sinking quality - local (same as global)+prop_quality_dec_l mx g (E (Seq h d _)) = let +    q10 = Just $ B.map (const 10) d+    q20 = Just $ B.map (const 20) d+    in Q.local_score mx g (Seq h d q10) (Seq h d q10) +           <= Q.local_score mx g (Seq h d q20) (Seq h d q20)++-- Sinking score with sinking quality - global+-- NB: score is always positive (aligning seq. with itself)+prop_quality_dec_g mx g (E (Seq h d _)) = let +    q10 = Just $ B.map (const 10) d+    q20 = Just $ B.map (const 20) d+    in Q.global_score mx g (Seq h d q10) (Seq h d q10)+           <= Q.global_score mx g (Seq h d q20) (Seq h d q20)++-- Sinking avg score with sinking quality +-- This is not correct, low quality reduces the penalty for mismatch+-- much harder than the reward for a match.  +-- (E.g. mmmxmmm may score higher with poor quality)+prop_quality_dec2 mx g (E (Seq h d _)) (E (Seq h' d' _)) = let +    q10 = Just $ B.map (const 10) d+    q20 = Just $ B.map (const 20) d+    (s1,a1) = Q.local_align mx g (Seq h d q10) (Seq h' d' q10) +    (s2,a2) = Q.local_align mx g (Seq h d q20) (Seq h' d' q20)+    in if null a1 && null a2 then True +       else s1 / fromIntegral (length a1) <= s2 / fromIntegral (length a2)++-- not expected to pass(?), but useful to debug indirect alignments+prop_indirect (E s1) (E s2) (E s3) = let a1 = snd $ A.global_align amx g s1 s2+                                         a2 = snd $ A.global_align amx g s2 s3+                                         a3 = snd $ A.global_align amx g s1 s3+                                         i1 = indirect a1 a2+                                     in if a3 == i1 then True +                                        else error ("\n"++unlines [showalign a1, showalign a2,"",showalign a3,"",showalign i1])++prop_indirect_recover (E s1) (E s2) (E s3) = +    let a1 = snd $ A.global_align amx g s1 s2+        a2 = snd $ A.global_align amx g s2 s3+        i1 = indirect a1 a2+        (s1',s3') = toStrings i1+        f = filter (/= '-')+    in if f s1' == toStr (seqdata s1) && f s3' == toStr (seqdata s3)+       then True else error ("\n"++(unlines [showalign a1, showalign a2,"",showalign i1]))++{-+-- | Affine = Simple when go == ge+prop_affine1 s1 s2 = S.global_align mx g s1 s2 == A.global_align mx (g,g) s1 s2+prop_affine2 s1 s2 = S.global_score mx g s1 s2 == A.global_score mx (g,g) s1 s2+-}++{-+-- | For a given scoring scheme, scoring must correspond to alignments+prop_score1 mx g (Eq s1) (Eq s2) = +    sum . S.eval mx g . snd . S.global_align mx g s1 $ s2 +    == S.global_score mx g s1 s2+-- -> todo: ditto for A and Q, ditto for local+-}++{-+-- | Any global alignment should return the exact same sequences+prop_recover mx g s1 s2 = undefined++-- | Qual = Affine when no quality data are supplied+prop_qual s1 s2 = undefined++prop_local s1 s2 = S.global_score mx g s1 s2 <= S.local_score mx g s1 s2+-- -> todo: ditto for A and Q++-- | Verify that all columns have the same lenght (ie. it should be a square matrix)+--   (Check vs sequence lengths)  Both QAlign.columns and AlignData.columns+prop_columns = undefined+-}++instance Arbitrary Soap.SoapAlign where+    arbitrary = elements [18..36] >>= genSoapAlign++genSeqName :: Gen BC.ByteString+genSeqName = liftM (BC.pack . filter (not . isSpace)) $ sized $ vector . (+ 1)++genSoapAlign :: Offset -> Gen Soap.SoapAlign+genSoapAlign len = do name <- genSeqName+                      sequ <- liftM (fromStr . map fromN) $ vector $ fromIntegral len +                      qual <- liftM (B.pack . map fromQ) $ vector $ fromIntegral len +                      let nhit = 1+                      let pairend = 'a'+                      strand <- elements [Fwd, RevCompl]+                      refname <- genSeqName+                      refstart <- genOffset+                      nmismatch <- choose (0, 3)+                      mismatches <- genMismatches nmismatch sequ qual+                      return $ Soap.SA name sequ qual nhit pairend len strand refname refstart nmismatch mismatches++genMismatches :: Int -> SeqData -> QualData -> Gen [Soap.SoapAlignMismatch]+genMismatches nmismatch sequ qual = shuffle [0..(BC.length sequ - 1)] >>= mapM genMismatchAt . take nmismatch+    where shuffle :: [a] -> Gen [a]+          shuffle = liftM (map fst . sortBy (comparing snd)) . mapM prioritize+          prioritize :: a -> Gen (a, Double)+          prioritize x = liftM ((,) x) $ choose (0.0, 1.0)+          genMismatchAt :: Offset -> Gen Soap.SoapAlignMismatch+          genMismatchAt off = let readnt = BC.index sequ off+                                  qualnt = B.index qual off+                              in do refnt <- elements $ filter (/= readnt) "acgt"+                                    return $ Soap.SAM readnt refnt off qualnt++prop_SoapAlign_invertparse :: Soap.SoapAlign -> Bool+prop_SoapAlign_invertparse sa = let saLine = Soap.unparse sa+                                    saErr :: Either String Soap.SoapAlign+                                    saErr = Soap.parse saLine+                                in case saErr of+                                     (Right sa') -> let saLine' = Soap.unparse sa'+                                                    in (sa == sa') && (saLine == saLine')+                                     (Left x) -> error x++genNumberedEsts :: Gen EST_set+genNumberedEsts = liftM (ESet . zipWith numberedSequence [1..]) $ replicateM numEsts genEstSeqData+    where numberedSequence z s = Seq (BC.pack $ "seq" ++ show z) s Nothing+          numEsts = 100+          genEstSeqData = liftM BC.pack $ choose (200, 1000) >>= flip replicateM (elements "ACGT")++genReadCSeqLoc :: Offset -> Sequence t -> Gen SeqLoc.ContigSeqLoc+genReadCSeqLoc len (Seq seqName seqSeq _) +    = do strand <- elements [Fwd, RevCompl]+         offset5 <- liftM fromIntegral $ chooseInteger (0, fromIntegral $ BC.length seqSeq - len)+         return $ OnSeq seqName $ CLoc.ContigLoc offset5 len strand+    where chooseInteger :: (Integer, Integer) -> Gen Integer+          chooseInteger = choose++genReadSequence :: Offset -> Sequence t -> Gen (Sequence t)+genReadSequence len sequ@(Seq _ seqSeq _) = do rloc <- genReadCSeqLoc len sequ+                                               return $ Seq (readName rloc) (readSeq rloc) (readQual rloc)+    where readName (OnSeq targName (CLoc.ContigLoc targOffset5 targLen targStrand))+              = BC.intercalate (BC.singleton '_') [ targName+                                                  , BC.pack $ show targOffset5+                                                  , BC.pack $ show targLen+                                                  , BC.pack $ show targStrand ]+          readSeq (OnSeq _ cloc) = either error id $ CLoc.seqData seqSeq cloc+          readQual rloc = Just $ B.replicate (BC.length $ readSeq rloc) (fromIntegral 40)++genEstReadSequence :: Offset -> EST_set -> Gen (Sequence Nuc)+genEstReadSequence len (ESet seqs) = elements seqs >>= genReadSequence len++genEstReads :: Offset -> EST_set -> Gen [Sequence Nuc]+genEstReads len eset = choose readNumRange >>= flip replicateM (genEstReadSequence len eset)+    where readNumRange = (5, 50)++runSoap :: EST_set -> [Sequence Nuc] -> IO [Soap.SoapAlign]+runSoap (ESet refSeqs) readSeqs+    = withTempFile "/tmp" "soap-refs-XXXXXX.txt"   $ \(refFile, hRef) ->+      withTempFile "/tmp" "soap-reads-XXXXXX.txt"  $ \(readFile, hRead) ->+      withTempFile "/tmp" "soap-aligns-XXXXXX.txt" $ \(alignFile, hAlign) ->+      do hWriteFasta hRef refSeqs >> hClose hRef+         hWriteFastQ hRead readSeqs >> hClose hRead+         hClose hAlign+         rawSystem soapExe [ "-d", refFile, "-a", readFile, "-r", "2", "-v", "3", "-o", alignFile ]+         liftM (map parseSoapLine . BC.lines) $ BC.readFile alignFile+    where parseSoapLine = either error id . Soap.parse+          soapExe = "/home/ingolia/Prog/extsrc/soap_1.10/soap.contig"+          withTempFile :: FilePath -> String -> ((FilePath,Handle) -> IO r) -> IO r+          withTempFile name tmp = bracket (openTempFile name tmp)+                                  (\(f,h) -> hClose h >> removeFile f)+++verifySoap :: [Soap.SoapAlign] -> Sequence t -> Bool+verifySoap aligns (Seq readName _ _)+    = case BC.split '_' readName of+        [refName, offset5Str, lenStr, strandStr]+            -> let offset5 = read $ BC.unpack offset5Str+                   len = read $ BC.unpack lenStr+                   strand = read $ BC.unpack strandStr+                   cLoc = CLoc.ContigLoc offset5 len strand+                   readLoc = OnSeq refName cLoc+               in any ((== readLoc) . Soap.refCSeqLoc) aligns+        _ -> False++prop_soap_in_vivo :: Property+prop_soap_in_vivo = forAll genNumberedEsts $ \refSeqs ->+                    forAll (genEstReads 36 refSeqs) $ \readSeqs ->+                    let aligns = unsafePerformIO $ runSoap refSeqs readSeqs+                    in collect (length aligns - length readSeqs) $ all (verifySoap aligns) readSeqs++bench :: [Test]+--           .........o.........o.........o+bench = [ T "local score, short" bench_local_rev_s+        , T "overlap score, short" bench_overlap_rev_s+        , T "global score, short" bench_global_rev_s++        , T "local score, long" bench_local_rev_l+        , T "overlap score, long" bench_overlap_rev_l+        , T "global score, long" bench_global_rev_l+        ]++bench_global_rev_s (ES x) (ES y) = +    abs (Q.global_score mx g x y+         - Q.global_score mx g (revcompl x) (revcompl y)) < 0.1++bench_local_rev_s (ES x) (ES y) = +    abs (Q.local_score mx g x y+         - Q.local_score mx g (revcompl x) (revcompl y)) < 0.1++bench_overlap_rev_s (ES x'@(Seq h d _)) (ES y'@(Seq h2 d2 _)) = +    let x = Seq h d Nothing+        y = Seq h2 d2 Nothing+        f = Q.overlap_align mx g x y+        r = Q.overlap_align mx g (revcompl x) (revcompl y)+    in if abs (fst f-fst r) < 0.1 then True+    else error ("\nFwd: "++show (fst f)++"\n"++ showalign (snd f)+                   ++"\nRev: "++show (fst r)++"\n"++ showalign (snd r)+                   ++"\n"++toStr (seqdata x)++"\n"++(toStr $ seqdata y))++bench_overlap_rev_l (EL x) (EL y) = +    abs (Q.overlap_score mx g x y+         - Q.overlap_score mx g (revcompl x) (revcompl y)) < 0.1++bench_global_rev_l (EL x) (EL y) = +    abs (Q.global_score mx g x y+         - Q.global_score mx g (revcompl x) (revcompl y)) < 0.1++bench_local_rev_l (EL x) (EL y) = +    abs (Q.local_score mx g x y+         - Q.local_score mx g (revcompl x) (revcompl y)) < 0.1

@@ -0,0 +1,38 @@+-- | Clustering Tests+module Bio.Clustering.Test where++import Bio.Util.TestBase+import Bio.Clustering+import Data.List (nub,sort)++tests :: [Test]+--           .........o.........o.........o+tests = [ T "retains elements"              prop_retains+        , T "hierarchy w/sorted"            prop_hierarchy+        , T "triangle ineq"                 prop_triangle+        ]++-- | Check that all elements from pairs are in the clustering+prop_retains :: [(Double,Int,Int)] -> Bool+prop_retains xs = clusterElements (cluster_sl xs) == listElements xs+    where listElements = nub . sort . concatMap (\(_,x,y)->[x,y])+          clusterElements = nub . sort . concatMap cE+          cE (Branch _ left right) = cE left ++ cE right+          cE (Leaf a) = [a]++-- | Check that the order of branches is correct, as long as the order of+--   input pairs are sorted.          +prop_hierarchy :: [Int] -> [Int] -> Bool+prop_hierarchy xs ys = let ts = zip3 [(1::Double)..] xs ys+                           cs = cluster_sl ts +                           isSorted (Leaf _) = True+                           isSorted (Branch s left right) = +                               lessThan s left && lessThan s right &&+                                        isSorted left && isSorted right+                           lessThan x (Branch y _ _) = x >= y+                           lessThan _ (Leaf _) = True+                       in all isSorted cs++prop_triangle :: [(Double,Int,Int)] -> [(Double,Int,Int)] -> Bool+prop_triangle xs ys = length (cluster_sl (xs ++ ys))+                      <= length (cluster_sl xs) + length (cluster_sl ys)

@@ -0,0 +1,64 @@+module Bio.GFF3.Test (tests)+    where++import Control.Monad+import Control.Monad.Error+import qualified Data.ByteString.Lazy.Char8 as LBS+import Data.Char+import Test.QuickCheck++import Bio.GFF3.Escape+import qualified Bio.GFF3.Feature as F+import Bio.Location.Strand+import Bio.Sequence.SeqData+import Bio.Util.TestBase++tests :: [Test]+tests = [ T "%-escape inversion"            test_Escape_invert+        , T "%-escape completeness"         test_Escape_completeness++        , T "GFF3 parse inversion"          test_GFF3_unparseParse+        ]++instance Arbitrary LBS.ByteString where+    arbitrary = liftM LBS.pack $ arbitrary++-- Bio.GFF3.Escape++test_Escape_invert :: [Char] -> LBS.ByteString -> Bool+test_Escape_invert escchrs str = let escstr = escapeAllOf ('%':escchrs) str+                                     unescstr :: Either String LBS.ByteString+                                     unescstr = unEscapeByteString escstr+                                 in unescstr == Right str++test_Escape_completeness :: [Char] -> LBS.ByteString -> Bool+test_Escape_completeness escchrs str+    = let allescchrs = ('%':escchrs)+          nescchrs = length $ LBS.findIndices (flip elem allescchrs) str          +      in (length $ LBS.elemIndices '%' $ escapeAllOf allescchrs str) == nescchrs++-- Bio.GFF3.Feature++instance Arbitrary F.GFFAttr where+    arbitrary = liftM2 F.GFFAttr arbitrary (sized vector)++instance Arbitrary Strand where+    arbitrary = elements [Fwd, RevCompl]++instance Arbitrary F.Feature where+    arbitrary = do seqid <- arbitrary+                   source <- arbitrary+                   stype <- arbitrary+                   seq5 <- genNonNegOffset+                   len <- genPositiveOffset+                   score <- return Nothing -- Tricky to ask strict equality in a Double+                   strand <- arbitrary+                   phase <- oneof [return Nothing, liftM (Just . fromIntegral) $ choose ((0::Int),2)]+                   attrs <- sized vector+                   return $ F.Feature seqid source stype seq5 (seq5 + len - 1) score strand phase attrs++test_GFF3_unparseParse :: F.Feature -> Bool+test_GFF3_unparseParse f = let fline = F.unparse f+                               f' :: Either String F.Feature+                               f' = F.parse fline+                           in f' == Right f

@@ -0,0 +1,326 @@+module Bio.Location.Test (tests)+    where++import Control.Monad+import Control.Monad.Error+import qualified Data.ByteString.Lazy.Char8 as LBS+import Data.Char+import Data.Either+import Data.Int (Int64)+import Data.Ix (inRange)+import Data.List+import Data.Maybe+import Test.QuickCheck++import qualified Bio.Location.ContigLocation as CLoc+import qualified Bio.Location.LocMap as LM+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 qualified Bio.Location.SeqLocMap as SLM+import Bio.Location.Strand+import Bio.Sequence.SeqData+import Bio.Util.TestBase++++tests :: [Test]+tests = [ T "Strand revCompl"               test_Strand_revCompl +        , T "Char revCompl"                 test_Char_revCompl+        , T "SeqData revCompl"              property_SeqData_revCompl+        , T "Sequence revCompl"             property_Sequence_revCompl++        , T "Pos revCompl"                  test_Pos_revCompl+        , T "Pos seqNt"                     property_Pos_seqNt+        , T "Pos seqNtPadded"               property_Pos_seqNtPadded++        , T "Contig revCompl"               test_Contig_RevCompl+        , T "Contig into/outof inversion"   property_ContigIntoOutof+        , T "Contig outof/into inversion"   property_ContigOutofInto+        , T "Contig into based on bounds"   test_Contig_IntoBounds+        , T "Contig outof based on bounds"  test_Contig_OutofBounds+        , T "Contig seqData"                property_Contig_seqData+        , T "Contig seqDataPadded"          property_Contig_seqDataPadded+        , T "Contig extend/revCompl"        property_Contig_extendRevCompl+        , T "Contig fromStartEnd"           property_Contig_fromStartEnd++        , T "Loc revCompl"                  test_Loc_RevCompl+        , T "Loc into/outof inversion"      property_LocIntoOutof+        , T "Loc outof/into inversion"      property_LocOutofInto+        , T "Loc outof based on bounds"     test_Loc_OutofBounds+        , T "Loc within"                    property_Loc_Within+        , T "Loc seqData"                   property_Loc_seqData+        , T "Loc seqDataPadded"             property_Loc_seqDataPadded++        , T "LocMap Within"                 property_LocMap_Within+        , T "LocMap Overlaps"               property_LocMap_Overlaps++        , T "SeqLocMap Within"              property_SeqLocMap_Within+        , T "SeqLocMap Overlaps"            property_SeqLocMap_Overlaps++        ]+++-- Bio.Location.Stranded++genNtSeqData :: Int -> Gen SeqData+genNtSeqData = liftM LBS.pack . flip replicateM (elements "ACGT")++test_revCompl :: (Eq s, Stranded s) => s -> Bool+test_revCompl s = (revCompl . revCompl) s == s++test_Strand_revCompl :: Strand -> Bool+test_Strand_revCompl = test_revCompl++test_Char_revCompl :: Char -> Bool+test_Char_revCompl = test_revCompl++property_SeqData_revCompl :: Property+property_SeqData_revCompl = forAll (sized genNtSeqData) test_revCompl++property_Sequence_revCompl :: Property+property_Sequence_revCompl+    = forAll arbitrary $ \name ->+      let mkSeq s = Seq name s Nothing+      in forAll (sized genNtSeqData) $ \sequ ->+          ((revcompl . mkSeq) sequ) == ((mkSeq . revCompl) sequ)++-- Bio.Location.Position++test_Pos_revCompl :: Pos.Pos -> Bool+test_Pos_revCompl = test_revCompl++property_Pos_seqNt :: Pos.Pos -> Property+property_Pos_seqNt pos@(Pos.Pos off str)+    = let pos = Pos.Pos off str+      in forAll genPositiveOffset $ \seqlen ->+         forAll (genNtSeqData $ fromIntegral seqlen) $ \sequ ->+         let actual = Pos.seqNt sequ pos+         in if inRange (0, seqlen - 1) off +            then let fwdNt = LBS.index sequ off+                 in if str == Fwd +                    then actual == Right fwdNt+                    else actual == Right (compl fwdNt)+            else isLeft actual+    where isLeft :: Either String Char -> Bool+          isLeft = either (const True) (const False)++property_Pos_seqNtPadded :: Pos.Pos -> Property+property_Pos_seqNtPadded pos@(Pos.Pos off str)+    = forAll genPositiveOffset $ \seqlen ->+      forAll (genNtSeqData $ fromIntegral seqlen) $ \sequ ->+          (Pos.seqNt sequ pos `catchError` returnN) == (Right $ Pos.seqNtPadded sequ pos)+    where returnN :: String -> Either String Char+          returnN _ = return 'N'++-- Bio.Location.ContigLocation++instance Arbitrary Strand where+    arbitrary = elements [Fwd, RevCompl]++instance Arbitrary Pos.Pos where+    arbitrary = liftM2 Pos.Pos genOffset arbitrary++instance Arbitrary CLoc.ContigLoc where+    arbitrary = liftM3 CLoc.ContigLoc genOffset genPositiveOffset arbitrary++instance Arbitrary LBS.ByteString where+    arbitrary = liftM LBS.pack $ arbitrary++test_Contig_RevCompl :: CLoc.ContigLoc -> Bool+test_Contig_RevCompl = test_revCompl++property_ContigIntoOutof :: CLoc.ContigLoc -> Pos.Pos -> Property+property_ContigIntoOutof contig pos+    = let !mInpos = CLoc.posInto pos contig+          !mOutpos = mInpos >>= flip CLoc.posOutof contig                     +      in (isJust mInpos) ==> mOutpos == (Just pos)++property_ContigOutofInto :: CLoc.ContigLoc -> Pos.Pos -> Property+property_ContigOutofInto contig pos+    = let !mOutpos = CLoc.posOutof pos contig+          !mInpos = mOutpos >>= flip CLoc.posInto contig+      in (isJust mOutpos) ==> mInpos == (Just pos)++test_Contig_IntoBounds :: CLoc.ContigLoc -> Pos.Pos -> Bool+test_Contig_IntoBounds contig pos+    = let !mInpos = CLoc.posInto pos contig+          !offset = Pos.offset pos+          !(cstart, cend) = CLoc.bounds contig+      in (isJust mInpos) == (offset >= cstart && offset <= cend)++test_Contig_OutofBounds :: CLoc.ContigLoc -> Pos.Pos -> Bool+test_Contig_OutofBounds contig pos+    = let !offset = Pos.offset pos+      in (isJust $ CLoc.posOutof pos contig) == (offset >= 0 && offset < CLoc.length contig)++property_Contig_seqData :: CLoc.ContigLoc -> Property+property_Contig_seqData contig+    = forAll (genNonNegOffset >>= genNtSeqData . fromIntegral) $ \sequ ->+      let seqData :: Either String SeqData+          seqData = CLoc.seqData sequ contig+          padded = CLoc.seqDataPadded sequ contig+      in case seqData of+           (Right subsequ) -> and [ padded == subsequ, 'N' `LBS.notElem` padded ]+           (Left _) -> 'N' `LBS.elem` padded++property_Contig_seqDataPadded :: CLoc.ContigLoc -> Property+property_Contig_seqDataPadded contig+    = forAll (genNonNegOffset >>= genNtSeqData . fromIntegral) $ \sequ ->+      (LBS.pack $ map (Pos.seqNtPadded sequ) contigPoses) == CLoc.seqDataPadded sequ contig+    where contigPoses = mapMaybe (flip CLoc.posOutof contig . flip Pos.Pos Fwd) [0..(CLoc.length contig - 1)]++property_Contig_extendRevCompl :: CLoc.ContigLoc -> Property+property_Contig_extendRevCompl contig+    = forAll (liftM2 (,) genNonNegOffset genNonNegOffset) $ \(ext5, ext3) ->+      (revCompl $ CLoc.extend (ext5, ext3) contig) == (CLoc.extend (ext3, ext5) $ revCompl contig)++property_Contig_fromStartEnd :: CLoc.ContigLoc -> Property+property_Contig_fromStartEnd contig+    = (CLoc.length contig > 1) ==>+      (CLoc.fromStartEnd (Pos.offset $ CLoc.startPos contig) (Pos.offset $ CLoc.endPos contig)) == contig++-- Bio.Location.Location++genInvertibleLoc :: Gen Loc.Loc+genInvertibleLoc = sized $ \sz -> do ncontigs <- choose (1, sz + 1)+                                     fwdloc <- liftM Loc.Loc $ genContigs ncontigs+                                     rc <- arbitrary+                                     if rc then return $ revCompl fwdloc else return fwdloc+    where genContigs = liftM (reverse . foldl' intervalsToContigs []) . genIntervals+          genIntervals nints = replicateM nints $ liftM2 (,) genPositiveOffset genPositiveOffset+          intervalsToContigs [] (init5, len) = [CLoc.ContigLoc init5 len Fwd]+          intervalsToContigs prevs@(prev:_) (nextoffset, nextlen)+              = let !prevend = CLoc.offset5 prev + CLoc.length prev+                in (CLoc.ContigLoc (prevend + nextoffset) nextlen Fwd):prevs++instance Arbitrary Loc.Loc where+    arbitrary = sized $ \sz -> do nintervals <- choose (1, sz + 1)+                                  liftM Loc.Loc $ vector nintervals++test_Loc_RevCompl :: Loc.Loc -> Bool+test_Loc_RevCompl = test_revCompl++property_LocIntoOutof :: Loc.Loc -> Pos.Pos -> Property+property_LocIntoOutof loc pos+    = let !mInpos = Loc.posInto pos loc+          !mOutpos = mInpos >>= flip Loc.posOutof loc+      in (isJust mInpos) ==> mOutpos == (Just pos)++property_LocOutofInto :: Pos.Pos -> Property+property_LocOutofInto pos+    = forAll genInvertibleLoc $ \loc ->+      let !mOutpos = Loc.posOutof pos loc+          !mInpos = mOutpos >>= flip Loc.posInto loc+      in (isJust mOutpos) ==> mInpos == (Just pos)++test_Loc_OutofBounds :: Loc.Loc -> Pos.Pos -> Bool+test_Loc_OutofBounds loc pos+    = let !offset = Pos.offset pos+      in (isJust $ Loc.posOutof pos loc) == (offset >= 0 && offset < Loc.length loc)++property_Loc_seqData :: Loc.Loc -> Property+property_Loc_seqData loc+    = forAll (genNonNegOffset >>= genNtSeqData . fromIntegral) $ \sequ ->+      let seqData :: Either String SeqData+          seqData = Loc.seqData sequ loc+          padded = Loc.seqDataPadded sequ loc+      in case seqData of+           (Right subsequ) -> and [ padded == subsequ, 'N' `LBS.notElem` padded ]+           (Left _) -> 'N' `LBS.elem` padded++property_Loc_seqDataPadded :: Loc.Loc -> Property+property_Loc_seqDataPadded loc+    = forAll (genNonNegOffset >>= genNtSeqData . fromIntegral) $ \sequ ->+      (LBS.pack $ map (Pos.seqNtPadded sequ) locPoses) == Loc.seqDataPadded sequ loc+    where locPoses = mapMaybe (flip Loc.posOutof loc . flip Pos.Pos Fwd) [0..(Loc.length loc - 1)]+++property_Loc_Within :: Pos.Pos -> Property+property_Loc_Within pos+    = forAll genInvertibleLoc $ \loc ->+      (pos `Loc.isWithin` loc) == (maybe False ((/= RevCompl) . Pos.strand) $ Loc.posInto pos loc)++--++class Checkable a where+    addCheck :: a -> a++instance Checkable () where addCheck = id+instance Checkable Int where addCheck = id+instance Checkable Char where addCheck = id++instance (Checkable a, Checkable b) => Checkable (a, b) where +    addCheck (x, y) = (addCheck x, addCheck y)++instance (Checkable a) => Checkable [a] where+    addCheck = map addCheck ++instance (Checkable a, Checkable b) => Checkable (a -> b) where+    addCheck f = \x -> (addCheck . f) (addCheck x)++instance Checkable (LM.LocMap a) where+    addCheck x = case LM.checkInvariants x of+                   [] -> x+                   errs -> error $ unlines errs++instance Checkable Int64 where addCheck = id+instance Checkable Pos.Pos where addCheck = id+instance Checkable Loc.Loc where addCheck = id++genLocs :: Gen [Loc.Loc]+genLocs = sized $ \sz -> choose (0, sz) >>= vector++property_LocMap_Within :: Pos.Pos -> Property+property_LocMap_Within seqpos +    = forAll genLocs $ \locs ->+      forAll genPositiveOffset $ \zonesize ->+      let !locents = zip locs ['0'..]+          !locmap = (addCheck LM.fromList) zonesize locents+          !hits = filter (Loc.isWithin seqpos . fst) locents+          !maphits = (addCheck LM.lookupWithin) seqpos locmap+      in -- collect (length hits) $ +         sort hits == sort maphits++property_LocMap_Overlaps :: Loc.Loc -> Property+property_LocMap_Overlaps loc+    = forAll genLocs $ \locs ->+      forAll genPositiveOffset $ \zonesize ->+      let !locents = zip locs ['0'..]+          !locmap = (addCheck LM.fromList) zonesize locents+          !hits = filter (Loc.overlaps loc . fst) locents+          !maphits = (addCheck LM.lookupOverlaps) loc locmap+      in -- collect (length hits) $ +        sort hits == sort maphits++--++genSeqs :: Gen [SeqName]+genSeqs = sized $ \sz -> choose (1, sz + 1) >>= vector++genSeqLocs :: [SeqName] -> Gen [SeqLoc.SeqLoc]+genSeqLocs seqNames = genLocs >>= mapM genSeqLoc+    where genSeqLoc loc = liftM (flip OnSeq loc) $ elements seqNames++property_SeqLocMap_Within+    = forAll genSeqs $ \seqs -> +      forAll (genSeqLocs seqs) $ \slocs ->+      forAll (liftM2 OnSeq (elements seqs) arbitrary) $ \spos ->+      let !slocents = zip slocs ['0'..]+          !slocmap = SLM.fromList slocents+          !hits = filter (andSameSeq Loc.isWithin spos . fst) slocents+          !maphits = SLM.lookupWithin spos slocmap+      in -- collect (length hits, length slocents) $ +        sort hits == sort maphits++property_SeqLocMap_Overlaps+    = forAll genSeqs $ \seqs ->+      forAll (genSeqLocs seqs) $ \slocs ->+      forAll (liftM2 OnSeq (elements seqs) arbitrary) $ \sloc ->+      let !slocents = zip slocs ['0'..]+          !slocmap = SLM.fromList slocents+          !hits = filter (andSameSeq Loc.overlaps sloc . fst) slocents+          !maphits = SLM.lookupOverlaps sloc slocmap+      in -- collect (length hits, length slocents) $ +        sort hits == sort maphits

@@ -78,7 +78,7 @@ import System.Directory (doesFileExist)  -- | Read nucleotide sequences in any format - Fasta, SFF, FastQ, 2bit, PHD...---   Todo: read quality automatically if available+--   Todo: detect Illumina vs Sanger FastQ, transparent compression readNuc :: FilePath -> IO [Sequence Nuc] readNuc fp     | ext `elem` ["fasta", "fna", "fa", "fst"] = do @@ -89,7 +89,8 @@        return (map castSeq ss)   | ext == "2bit"                     = read2Bit                               $ fp   | ext == "sff"                      = fmap sffToSequence . readSFF           $ fp-  | ext `elem` ["fq","fastq"]         = readFastQ                              $ fp+  | ext `elem` ["fq","fastq"]         = readSangerQ                            $ fp+  | ext == "txt"                      = readIllumina                           $ fp   | ext2 == "phd"                     = fmap return . readPhd                  $ fp -- only a single sequence   -- "ace" ?   | otherwise                         = error "readNuc: unknown file suffix!"

@@ -28,7 +28,7 @@ ) where  -import Data.Char (chr) -- isSpace+import Data.Char (chr, isSpace) import Data.List (groupBy,intersperse) import System.IO @@ -59,12 +59,7 @@  -- | Lazily read sequences from a FASTA-formatted file readFasta :: FilePath -> IO [Sequence Unknown]-readFasta f = (mkSeqs . blines) `fmap` B.readFile f---- | Replacement for 'lines' that gobbles control-M's--- Some tools, like CLC, likes to add these to the end of each line.-blines :: B.ByteString -> [B.ByteString]-blines = B.lines . B.filter (/=Data.Char.chr 13)+readFasta f = (mkSeqs . B.lines) `fmap` B.readFile f  -- | Write sequences to a FASTA-formatted file. --   Line length is 60.@@ -76,7 +71,7 @@  -- | Read quality data for sequences to a file. readQual :: FilePath -> IO [Sequence Unknown]-readQual f = (mkQual . blines) `fmap` B.readFile f+readQual f = (mkQual . B.lines) `fmap` B.readFile f  -- | Write quality data for sequences to a file. writeQual :: FilePath -> [Sequence a] -> IO ()@@ -120,7 +115,7 @@  -- | Lazily read sequence from handle hReadFasta :: Handle -> IO [Sequence Unknown]-hReadFasta h = (mkSeqs . blines) `fmap` B.hGetContents h+hReadFasta h = (mkSeqs . B.lines) `fmap` B.hGetContents h  -- | Write sequences in FASTA format to a handle. hWriteFasta :: Handle -> [Sequence a] -> IO ()@@ -161,7 +156,7 @@ mkSeq (l:ls)    -- maybe check this?  | B.length l < 2 || isSpace (B.head $ B.tail l)    --  = error "Trying to read sequence without a name...and failing."-  | otherwise = Seq (B.drop 1 l) (B.concat $ takeWhile isSeq ls) Nothing+  | otherwise = Seq (B.drop 1 l) (B.filter (not . isSpace) $ 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" @@ -185,6 +180,6 @@ countSeqs :: FilePath -> IO Int countSeqs f = do   ss <- B.readFile f-  let hdrs = filter (('>'==).B.head) $ filter (not . B.null) $ blines ss+  let hdrs = filter (('>'==).B.head) $ filter (not . B.null) $ B.lines ss   return (length hdrs) 

@@ -1,12 +1,14 @@-{- | Read (and write?) the SFF file format used by+{- | Read and write the SFF file format used by    Roche\/454 sequencing to store flowgram data.     A flowgram is a series of values (intensities) representing homopolymer runs of    A,G,C, and T in a fixed cycle, and usually displayed as a histogram. -   The Staden Package contains an io_lib, with a C routine for parsing this format.+   This file is based on information in the Roche FLX manual.  Among other sources for information about+   the format, are The Staden Package, which contains an io_lib with a C routine for parsing this format.    According to comments in the sources, the io_lib implementation is based on a file-   called getsff.c, which I've been unable to track down.+   called getsff.c, which I've been unable to track down.  Other software parsing SFFs +   are QIIME, sff_extract, and Celera's sffToCa.     It is believed that all values are stored big endian. -}@@ -17,6 +19,7 @@                         , sffToSequence, rbToSequence                         , trim, trimFromTo, trimKey                         , baseToFlowPos, flowToBasePos+                        , trimFlows                         , test, convert, flowgram                         , masked_bases, cumulative_index                         , packFlows, unpackFlows@@ -49,67 +52,96 @@ type Index = Word8  -- Global variables holding static information+-- | An SFF file always start with this magic number. magic :: Int32 magic = 0x2e736666 +-- | Version is always 1. versions :: [Int32] versions = [1] +-- | Read an SFF file. readSFF :: FilePath -> IO SFF readSFF f = return . decode =<< LB.readFile f +-- | Extract the read without the initial (TCAG) key. 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)++")!") +-- | Extract the sequences from an 'SFF' data structure. sffToSequence :: SFF -> [Sequence Nuc] sffToSequence (SFF _ rs) = map rbToSequence rs +-- | Extract the sequence information from a 'ReadBlock'. rbToSequence :: ReadBlock -> Sequence Nuc rbToSequence r = clip (read_header r, bases r, quality r)     where 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)])+                           in 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) +-- | Trim a 'ReadBlock' limiting the number of flows.  If writing to+--   an SFF file, make sure you update the 'CommonHeader' accordingly.+--   See @examples/Flx.hs@ for how to use this.  +trimFlows :: Integral i => i -> ReadBlock -> ReadBlock+trimFlows l rb = rb { read_header = rh { num_bases = fromIntegral n+                                       , clip_qual_right = min cqr $ fromIntegral n+                                       }+                    , flow_data   = B.take (2*fromIntegral l) (flow_data rb)+                    , flow_index  = B.take n (flow_index rb)+                    , bases       = LB.take (fromIntegral n) (bases rb)+                    , quality     = LB.take (fromIntegral n) (quality rb)+                    }+  where n = (flowToBasePos rb l)-1+        rh = read_header rb+        cqr = clip_qual_right rh+ -- 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.++-- | Trim a read to specific sequence position, inclusive bounds+-- The current implementation has the unintended side effect of +-- always trimming the flowgram down to a basecalled position.  +-- Note that you can't (easily) write trimmed 'ReadBlock's to a file,+-- since they need to have the same number of flows as given in +-- the 'CommmonHeader'. 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)-                        rh = read_header rd-                        [r',l'] = map fromIntegral [r,l]-                        rh' = rh { num_bases = fromIntegral (r'-l'+1)-                                 , clip_qual_left = max 0 $ clip_qual_left rh-l'-                                 , clip_qual_right = min (clip_qual_right rh-l') (r'-l'+1)-                                 }-                    in rd { read_header = rh'-                          , flow_data = trim_flw (flow_data rd)-                          , flow_index = trim_seq' (flow_index rd)-                          , bases = trim_seq (bases rd)-                          , quality = trim_seq (quality rd)-                          }+trimFromTo (l+1) 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)+  rh = read_header rd+  [r',l'] = map fromIntegral [r,l]+  rh' = rh { num_bases = fromIntegral (r'-l')+           , clip_qual_left = max 0 $ clip_qual_left rh-l'+           , clip_qual_right = min (clip_qual_right rh-l') (r'-l'+1)+           }+  in rd { read_header = rh'+        , 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+trim rb = let rh = read_header rb in trimFromTo (clip_qual_left rh) (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+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 +-- | Read an SFF file, but be resilient against errors. recoverSFF :: FilePath -> IO SFF recoverSFF f = return . unRecovered . decode =<< LB.readFile f @@ -130,6 +162,7 @@   hClose h   return $ fromIntegral c +-- | Write 'ReadBlock's to a file handle. writeReads :: Handle -> Int -> [ReadBlock] -> IO Int32 writeReads _ _ [] = return 0 writeReads h i (r:rs) = do@@ -192,6 +225,7 @@       put hd       mapM_ (put . RBI (fromIntegral $ flow_length hd)) rds +-- | Helper function for decoding a 'ReadBlock'. {-# INLINE getRB #-} getRB :: CommonHeader -> ReadHeader -> Get ReadBlock getRB chead rh = do@@ -236,17 +270,18 @@  -- ---------------------------------------------------------- -- | SFF has a 31-byte common header---   Todo: remove items that are derivable (counters, magic, etc)---   cheader_lenght points to the first read header.---   Also, the format is open to having the index anywhere between reads,+--+--   The format is open to having the index anywhere between reads, --   we should really keep count and check for each read.  In practice, it --   seems to be places after the reads. --    --   The following two fields are considered part of the header, but as --   they are static, they are not part of the data structure---     magic   :: Word32   -- ^ 0x2e736666, i.e. the string ".sff"---     version :: Word32   -- ^ 0x00000001-+--+-- @        +--     magic   :: Word32   -- 0x2e736666, i.e. the string \".sff\"+--     version :: Word32   -- 0x00000001+-- @ data CommonHeader = CommonHeader {           index_offset                            :: Int64    -- ^ Points to a text(?) section         , index_length, num_reads                 :: Int32@@ -288,7 +323,7 @@            }  -- ---------------------------------------------------------- --- | Each Read has a fixed read header+-- | Each Read has a fixed read header, containing various information. data ReadHeader = ReadHeader {       name_length                           :: Int16     , num_bases                             :: Int32@@ -335,9 +370,11 @@     , quality                    :: ! QualData     } +-- | Helper function to access the flowgram flowgram :: ReadBlock -> [Flow] flowgram = unpackFlows . flow_data +-- | Extract the sequence with masked bases in lower case masked_bases :: ReadBlock -> SeqData masked_bases rb = let   l = fromIntegral $ clip_qual_left $ read_header rb@@ -347,6 +384,7 @@                 , LBC.take r (LBC.drop (l-1) s)                 , LBC.map toLower $ LBC.drop r s] +-- | Extract the index as absolute coordinates, not relative. cumulative_index :: ReadBlock -> [Int] cumulative_index = scanl1 (+) . map fromIntegral . B.unpack . flow_index @@ -399,6 +437,7 @@   buf <- getLazyByteString 20   decodeSaneH prefix buf   +-- | Decode a 'ReadHeader', verifying that the data make sense. decodeSaneH :: String -> LBC.ByteString -> Get ReadHeader decodeSaneH prefix buf = do   let PartialReadHeader rhl nl _nb _ql _qr _al _ar rn = decode buf

@@ -1,80 +1,108 @@--- | This implements a number of filters used in the Titanium pipeline+-- | This implements a number of filters used in the Titanium pipeline, +--   based on published documentation. module Bio.Sequence.SFF_filters where -import Bio.Sequence.SFF (ReadBlock(..), ReadHeader(..), flowToBasePos, flowgram)+import Bio.Sequence.SFF (ReadBlock(..), ReadHeader(..)+                        , flowToBasePos, flowgram, cumulative_index) import qualified Data.ByteString.Lazy as B+import qualified Data.ByteString as SB import qualified Data.ByteString.Lazy.Char8 as BL import Data.List (tails)+import Data.Char (toUpper)  -- Ti uses a set of filters, described in the (something) manual. -- (GS Run Processor Application, section 3.2.2++) --- ** Discarding filters **+-- ** 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+-- | This filter discards empty sequences.+discard_empty :: DiscardFilter+discard_empty rb = num_bases (read_header rb) >= 5 -filter_key rb = ("TCAG"==) $ take 4 $ BL.unpack $ bases rb+-- | Discard sequences that don't have the given key tag (typically TCAG) at the start+--   of the read.+discard_key :: String -> DiscardFilter+discard_key key rb = (map toUpper key==) $ take (length key) $ 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)+--   before the last postitive flow.  The percentage can be given as a parameter.+discard_dots :: Double -> DiscardFilter+discard_dots p rb = let dotcount = SB.length $ SB.filter (>3) $ flow_index rb+                    in fromIntegral dotcount / fromIntegral (BL.length $ bases rb) < p+                       && last (cumulative_index rb) >= 84  -- | 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_mixed :: DiscardFilter+discard_mixed rb = let fs = dropWhile (<50) . reverse . flowgram $ rb+                       fl = dlength fs+                   in and+                      [ (dlength (filter (>50) fs) / fl) < 0.7 -- 70% positive+                      , (dlength (filter (<45) fs) / fl) > 0.3 -- 30% noise+                      , (dlength (filter (>75) fs) / fl) > 0.3 -- 30% postivie+                      , (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+discard_length :: Int -> DiscardFilter+discard_length n rb = length (flowgram rb) >= n --- ** Trimming filters **+-- ** 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)+trim_sigint :: TrimFilter+trim_sigint rb = clipSeq rb (sigint rb) +-- n counts the "bad" flow values, m counts flow position 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+sigint rb = let bs = drop 1 $ 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 +                xs = dropWhile (\(_,_,f) -> f<=70) +                     $ dropWhile (\(n,m,_)->(1000*n) `div` m > (30::Int)) +                     $ reverse bs+            in case xs of []          -> error "no sequence left?"+                          ((_,m,_):_) -> flowToBasePos rb m +-- | 3.2.2.1.5 Primer filter +-- This looks for the B-adaptor at the end of the read.  The 454 implementation isn't very+-- effective at finding mutated adaptors.+trim_primer :: String -> TrimFilter+trim_primer s rb = clipSeq rb (find_primer s rb) --- 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.+find_primer :: String -> ReadBlock -> Int+find_primer s rb = go (num_bases (read_header rb) - 10)+  where go i | i <= 5    = fromIntegral (num_bases $ read_header rb)+             | match i   = fromIntegral i+             | otherwise = go (i-1)+        match j = s' `B.isPrefixOf` B.drop (fromIntegral j) (bases rb)+        s' = BL.pack $ map toUpper $ take 14 s --- | 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+-- 3.2.2.1.6 Trimback valley filter is ignored, we don't understand the description. -qual20 :: ReadBlock -> Int-qual20 rs = (fromIntegral $ num_bases $ read_header rs) -            - (length . takeWhile (<20) . map (avg . take 10) . tails . reverse . B.unpack $ quality rs)+-- | 3.2.2.1.7 Quality score trimming trims using a 10-base window until a Q20 average is found.+trim_qual20 :: Int -> TrimFilter+trim_qual20 w rs = clipSeq rs $ qual20 w rs --- ** Utility functions **+qual20 :: Int -> ReadBlock -> Int+qual20 w rs = (fromIntegral $ num_bases $ read_header rs)+              - (length . takeWhile (<20) . map (avg . take w) . 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+avg xs = sum (map fromIntegral xs) / dlength xs  -- | Translate a number of flows to position in sequence, and update clipping data accordingly clipFlows :: ReadBlock -> Int -> ReadBlock@@ -85,3 +113,55 @@ 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 }}++-- ** Data++-- Celera docs, at http://sourceforge.net/apps/mediawiki/wgs-assembler/index.php?title=SffToCA++-- These are used for mate-pair libraries, should be located around the middle of the read:++flx_linker = "GTTGGAACCGAAAGGGTTTGAATTCAAACCCTTTCGGTTCCAAC"  -- Celera+ti_linker  = "TCGTATAACTTCGTATAATGTATGCTATACGAAGTTATTACG"  -- 20K cod jump++-- ti_linker and this: "AGCATATTGAAGCATATTACATACGATATGCTTCAATAATGC"+-- from "GS FLX Titanium 3 kb Span Paired End Library Preparation Method Manual April 2009"+-- ftp://ftp.genome.ou.edu/pub/for_broe/titanium/++-- These are used at the end of RNA (cDNA) sequences, after the poly-A tail:++rna_adapter   = "ggcgggcgatgtctcgtctgagcgggctggcaaggc" -- cod transcripts?+rna_adapter2  = "ttcgcagtgagtgacaggctagtagctgagcgggctggcaaggc"  -- Cod_c.sff+rna_adapter3  = "gacggggcggatgtctcgtctgagcgggcgtggcaaggc"       -- COD1.sff++-- These are used at the end of DNA sequencing reads:++rapid_adapter = "agtcgtggaggcaaggcacacagggatagg"  -- sea louse reads, key GACT+ti_adapter_b  = "ctgagactgccaaggcacacagggggatagg"  -- sea bass and l.s.Ca+                 +{- Email from Markus Grohme:+>> Here are the sequences as filed by 454 in their patent application:+>>> AdaptorA+>> CTGAGACAGGGAGGGAACAGATGGGACACGCAGGGATGAGATGG+>>> AdaptorB+>> CTGAGACACGCAACAGGGGATAGGCAAGGCACACAGGGGATAGG+>>+>> However, looking through some earlier project data I had, I also+>> retrieved the following (by simply making a consensus of+>> sequences that did not match the target genome anymore):+>>> 5prime454adaptor???+>> GCCTCCCTCGCGCCATCAGATCGTAGGCACCTGAAA+>>> 3prime454adaptor???+>> GCCTTGCCAGCCCGCTCAGATTGATGGTGCCTACAG+>>+>> I currently know one linker sequence (454/Roche also calls it spacer+>> for GS20 and FLX paired-end sequencing:+>>> flxlinker+>> GTTGGAACCGAAAGGGTTTGAATTCAAACCCTTTCGGTTCCAAC+>> For Titanium data using standard Roche protocol, you need to screen+>> for two linker sequences:+>>> titlinker1+>> TCGTATAACTTCGTATAATGTATGCTATACGAAGTTATTACG+>>> titlinker2+>> CGTAATAACTTCGTATAGCATACATTATACGAAGTTATACGA++-}

@@ -14,6 +14,9 @@       --   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.+      --+      --   If you use overloaded strings (e.g., ghc -XOverloadedString), you can+      --   easily construct sequences from string literals.       Sequence(..), Offset, SeqData,        -- | Quality data is normally associated with nucleotide sequences@@ -27,6 +30,7 @@     , appendHeader, setHeader        -- * Converting to and from [Char]+      -- | It is probably better to use the 'IsString' class from 'Data.String' for this.     , fromStr, toStr        -- * Sequence utilities

@@ -0,0 +1,119 @@+-- Test for sequence functionality++module Bio.Sequence.Test where++import Test.QuickCheck+import System.IO.Unsafe+import Data.Maybe (isJust)++import Bio.Sequence.HashWord+import Bio.Util.TestBase+import Bio.Sequence++tests_io :: [Test]+--              .........o.........o.........o+tests_io = [ T "serializing"  prop_serialize+           , T "serialize qual" prop_serialize_qual+           , T "serialize fasta+qual" prop_serialize_fastaQual+           , T "serialize fasta+qual multi" prop_serialize_fastaQual2+           , T "serialize fastq" prop_serialize_fastq+           , T "serialize fastq multi" prop_serialize_fastq_multi+           ]++prop_serialize (E s) = +    let [s'] = unsafePerformIO (do writeFasta "/tmp/serialize_test" [s]+                                   readFasta "/tmp/serialize_test")+    in s == castToNuc s'++prop_serialize_qual (Eq s@(Seq h d q)) = +    let [(Seq h' d' q')] = unsafePerformIO +                           (do writeQual "/tmp/serialize_qual" [s]+                               readQual "/tmp/serialize_qual")+    in h == h' && q == q'++prop_serialize_fastaQual (Eq s) = +    let [s'] = unsafePerformIO +                           (do writeFastaQual "/tmp/serialize_fasta" "/tmp/serialize_qual" [s]+                               readFastaQual "/tmp/serialize_fasta" "/tmp/serialize_qual")+    in  s == castToNuc s'++prop_serialize_fastaQual2 :: [ESTq] -> Bool+prop_serialize_fastaQual2 es = +    let ests = map (\(Eq x) -> x) es+        ests' = unsafePerformIO +                (do writeFastaQual "/tmp/serialize_fasta" "/tmp/serialize_qual" ests+                    readFastaQual "/tmp/serialize_fasta" "/tmp/serialize_qual")+    in ests == map castToNuc ests'++prop_serialize_fastq (Eq s) = +    let [s'] = unsafePerformIO +                           (do writeFastQ "/tmp/serialize_fastq" [s]+                               readFastQ "/tmp/serialize_fastq")+    in s' == s++prop_serialize_fastq_multi :: [ESTq] -> Bool+prop_serialize_fastq_multi es = +    let ests = map (\(Eq x) -> x) es+        ests' = unsafePerformIO +                (do writeFastQ "/tmp/serialize_fastq" ests+                    readFastQ "/tmp/serialize_fastq")+    in ests' == ests++-- ----------------------------------------------------------+-- Tests for HashWord++tests_hw :: [Test]++--              .........o.........o.........o+tests_hw = [ T "n2k vs k2n"   prop_n2k_k2n+           , T "contigous_0"  prop_contigous_0+           , T "prop_rcontig_0" prop_rcontig_0+           , T "prop_rcontig_1" prop_rcontig_1+           , T "prop_rclast" prop_rclast+           ]++prop_n2k_k2n :: Int -> Bool+prop_n2k_k2n i' = let i = abs i' `mod` 65536 in (n2k 8 . k2n 8) i == i++-- check that hashes is equal to hash over all indices+prop_contigous_0 k (E s) = k > 0 ==> hashes (contigous k) (seqdata s) == let indices =  [0..seqlength s-fromIntegral k] +                       in map (\(Just i,j)->(i,j)) $ filter (isJust.fst) $ zipWith (,) (map (hash (contigous k) (seqdata s)) indices) indices++-- rcontig is the minimum of hash of each forward word and its reverse complement+prop_rcontig_0 k (E s) = k > 0 ==> zipWith min (map fst . hashes (contigous k) . seqdata $ s) (map fst . reverse . hashes (contigous k) . seqdata . revcompl $ s)+                           == (map fst $ hashes (rcontig k) (seqdata s))++-- check that reverse (hashes . reverse) == id+prop_rcontig_1  k (E s) = k > 0 ==> (reverse . map fst . hashes (rcontig k) . seqdata . revcompl $ s) +                            == (map fst . hashes (rcontig k) . seqdata $ s)++-- remove duplicates, and check key values vs rcontig+prop_rcpacked_1 = undefined++-- last hash is equal to first hash on revcompl seq.+-- see hashcount below+prop_rclast k (E s) = k > 0 && (not . null . hs $ s) ==> rcl k s+    where hs = map fst . hashes (rcontig k) . seqdata++-- really only Nuc+rcl :: Int -> Sequence Nuc -> Bool+rcl k s = ((last . hs $ s) == (head . hs . revcompl $ s))+    where hs = map fst . hashes (rcontig k) . seqdata++-- benchmarks: todo: time hash generation for contig, rcontig, and gapped (when implemented)+bench = [ T "rc hash counts int  (8)" (hashcount_int  8)+        , T "rc hash counts int (16)" (hashcount_int 16)+        , T "rc hash counts     (16)" (hashcount 16)+        , T "rc hash counts     (32)" (hashcount 32)+        ]++hashcount, hashcount_int :: Int -> EST_set -> Property++hashcount k es' = k > 0 ==> let ESet es = es' +                                hs :: Sequence Nuc -> [Integer]+                                hs = map fst . hashes (rcontig k) . seqdata+                            in and $ map (\e -> null (hs e) || rcl k e || error (show k ++"\n" ++ toStr (seqdata e))) es+hashcount_int k es' = k > 0 ==> let ESet es = es' +                                    hs :: Sequence Nuc -> [Int]+                                    hs = map fst . hashes (rcontig k) . seqdata +                                in and $ map (\e -> null (hs e) || rcl k e || error (toStr $ seqdata e)) es

@@ -35,7 +35,7 @@ import Data.Bits  -- import Test.QuickCheck hiding (check)    -- QC 1.0-import Test.QuickCheck hiding ((.&.)) -- QC 2.0+-- import Test.QuickCheck hiding ((.&.)) -- QC 2.0  -- constants default_magic, default_version :: Word32

@@ -0,0 +1,19 @@+module Bio.Util.Test where++import Bio.Util as U+import Bio.Util.TestBase++import qualified Data.ByteString.Lazy.Char8 as BC++tests :: [Test]+tests = [ T "BS lines replacement" prop_bs_lines+        , T "BS internal lines"    prop_bs_mylines+        ]++-- test the 'lines' that is going to be used+prop_bs_lines :: String -> Bool+prop_bs_lines xs = Prelude.lines xs == (map BC.unpack . U.lines . BC.pack) xs++-- test mylines, which may be used or not, depending on whether we detect the LBS bug+prop_bs_mylines :: String -> Bool+prop_bs_mylines xs = Prelude.lines xs == (map BC.unpack . mylines . BC.pack) xs

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

@@ -0,0 +1,61 @@+GHCMAKE    = ghc --make -O+RUNHASKELL = runhaskell++default: +	@echo "This Makefile remains for historical reasons, and"+	@echo "For most purposes, you should use 'cabal-install' instead."+	@echo "But if you insist, please provide a make target"++clean:+	find . -name *.o -o -name *.hi -o -name *.p_o -o -name *.p_hi | xargs rm++user_install: user_conf build haddock user_inst++install: conf build haddock inst++user_conf:+	$(RUNHASKELL) Setup.hs configure -p --user --prefix=${HOME}++build:+	$(RUNHASKELL) Setup.hs build ++haddock: +	$(RUNHASKELL) Setup.hs haddock++user_inst:+	$(RUNHASKELL) Setup.hs install --user++conf:+	$(RUNHASKELL) Setup.hs configure -p ++inst:+	sudo $(RUNHASKELL) ./Setup.hs install++test:+	$(GHCMAKE) Test.hs -o qc+	./qc++test_hpc:+	$(GHCMAKE) -fhpc Test.hs -o qc.hpc+	./qc.hpc+	hpc report qc.hpc++# bench:+# 	$(GHCMAKE) Bench.hs -o qb+# 	@echo ==================== >> benchmark.log+# 	@echo -n Start: >> benchmark.log+# 	@date >> benchmark.log+# 	@uname -a >> benchmark.log+# 	@echo ghc version: `strings qb | grep '^[6-9]\.[0-9]'` >> benchmark.log+# 	@echo $(GHCMAKE) >> benchmark.log+# 	./qb +RTS -sbenchmark.gc | tee -a benchmark.log+# 	@echo -n End: >> benchmark.log+# 	@date >> benchmark.log++# bench_hpc:+# 	$(GHCMAKE) Bench.hs -o qb.hpc+# 	./qb.hpc+# 	hpc report qb.hpc++update:+	darcs pull http://malde.org/~ketil/biohaskell/biolib

@@ -0,0 +1,38 @@+-- Tests for all bio functionality++module Main where++import Test.QuickCheck+import System.IO++import Bio.Util.TestBase+import Bio.Util.Test as U+import Bio.Alignment.Test as A+import Bio.GFF3.Test as G+import Bio.Location.Test as L+import Bio.Sequence.Test as S+import Bio.Clustering.Test as C++all_tests :: [(String,[Test])]+all_tests = [ ("Sequence tests", S.tests_io)+            , ("Sequence hash tests", S.tests_hw)+            , ("Alignment tests", A.tests)+            , ("GFF3 tests", G.tests)+            , ("Location tests", L.tests)+            , ("Clustering tests", C.tests)+            , ("Util functions", U.tests)+            ]++main = do+  hSetBuffering stdout NoBuffering+  mapM_ runTests all_tests++runTests :: (String,[Test]) -> IO ()+runTests (hd,ts) = do+  putStrLn ("\n --- "++hd++" ---\n")+  mapM_ runTest ts++runTest (T name test) = do +  putStr $ name ++ " " ++ replicate (28-length name) '.' ++ " "+  quickCheck test+

@@ -1,7 +1,8 @@ Name:                bio-Version:             0.4.8+Version:             0.5 License:             LGPL License-file:        LICENSE+Cabal-Version:       >= 1.6 Author:              Ketil Malde Maintainer:          ketil@malde.org @@ -13,34 +14,43 @@                      .                      Current list of features includes: a Sequence data type supporting                      protein and nucleotide sequences and conversion between them.  As of version-		     0.4, different kinds of sequence have different types.  Support for quality+                     0.4, different kinds of sequence have different types.  Support for quality                      data, reading and writing Fasta formatted files, reading TwoBit 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.-		     .+                     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.+                     .                      The Darcs repository is at: <http://malde.org/~ketil/biohaskell/biolib>.-Homepage:            http://blog.malde.org/index.php/the-haskell-bioinformatics-library/+Homepage:            http://biohaskell.org/Libraries/Bio  Tested-With:         GHC==6.12.1 Build-Type:          Simple-Build-Depends:       base>=4 && <5, QuickCheck>=2, binary >=0.4 && <0.5, tagsoup>=0.8, bytestring >= 0.9.1,-                     containers, array, parallel, parsec, random, old-time, mtl, directory+Data-Files:          README, Makefile -Data-Files:          README+Flag test+  Description:       Build testing software+  Default:           False -Exposed-modules:     Bio.Sequence,+Flag examples+  Description:       Build example programs+  Default:           True++Library+  Build-Depends:     base>=4 && <5, binary >=0.4 && <0.5, tagsoup>=0.8, bytestring >= 0.9.1,+                     containers, array, parallel, parsec, mtl, directory++  Exposed-modules:   Bio.Sequence,                      Bio.Sequence.SeqData,                      Bio.Sequence.Fasta, Bio.Sequence.FastQ,-		     Bio.Sequence.TwoBit, Bio.Sequence.Phd,+                     Bio.Sequence.TwoBit, Bio.Sequence.Phd,                      Bio.Sequence.Entropy, Bio.Sequence.HashWord,                      Bio.Sequence.GOA,                      Bio.Sequence.GeneOntology,-		     Bio.Sequence.KEGG,-		     Bio.Sequence.AminoProperties,-		     Bio.Sequence.SFF, Bio.Sequence.SFF_name, Bio.Sequence.SFF_filters+                     Bio.Sequence.KEGG,+                     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,@@ -48,14 +58,76 @@                      Bio.Alignment.Multiple, Bio.Alignment.ACE,                      Bio.Alignment.Bowtie,                      Bio.Alignment.Soap,-		     Bio.Alignment.BED, Bio.Alignment.PSL+                     Bio.Alignment.BED, Bio.Alignment.PSL                      Bio.Clustering,                      Bio.Util, Bio.Util.Parsex-		 Bio.Location.Strand, Bio.Location.Position,-		 Bio.Location.ContigLocation, Bio.Location.Location, Bio.Location.LocMap,-		 Bio.Location.OnSeq, Bio.Location.SeqLocation, Bio.Location.SeqLocMap,-		 Bio.GFF3.Escape, Bio.GFF3.Feature, Bio.GFF3.FeatureHier, Bio.GFF3.FeatureHierSequences,-		 Bio.GFF3.SGD+                     Bio.Location.Strand, Bio.Location.Position,+                     Bio.Location.ContigLocation, Bio.Location.Location, Bio.Location.LocMap,+                     Bio.Location.OnSeq, Bio.Location.SeqLocation, Bio.Location.SeqLocMap,+                     Bio.GFF3.Escape, Bio.GFF3.Feature, Bio.GFF3.FeatureHier, Bio.GFF3.FeatureHierSequences,+                     Bio.GFF3.SGD -extensions:          CPP, ParallelListComp-ghc-options:         -Wall -O2 -fexcess-precision -funbox-strict-fields -auto-all+  Extensions:        CPP, ParallelListComp+  Ghc-Options:       -Wall -O2 -fexcess-precision -funbox-strict-fields -auto-all -fno-warn-unused-do-bind++Source-Repository head+  Type:           darcs+  Location:       http://malde.org/~ketil/biohaskell/biolib++Executable qc+  Main-Is:        Test.hs+  Other-Modules:  Bio.Util.TestBase, Bio.Util.Test,+                  Bio.Alignment.Test, Bio.GFF3.Test, Bio.Location.Test, Bio.Sequence.Test, Bio.Clustering.Test+  Hs-Source-Dirs: .  +  Build-Depends:  base >= 3 && <5, process, containers, random, QuickCheck >= 2, old-time+  if flag(test)+     Buildable:   True+  else+     Buildable:   False++-- Test.QuickBench needs to be ported to QC2.  Or everything moved over to criterion.+-- Executable qb+--  Main-Is:        Bench.hs+--  Other-Modules:  Test.QuickBench, Bio.Util.TestBase, Bio.Util.Test, +--                  Bio.Alignment.Test, Bio.Sequence.Test, Bio.Clustering.Test+--  Hs-Source-Dirs: .  +--  Build-Depends:  base >= 3 && <5, process, containers, random++Executable flx+  Main-Is:        Flx.hs+  Hs-Source-Dirs: examples .+  Other-Modules:  Bio.Sequence.SFF+  Build-Depends:  base >= 3 && <5, bytestring+  if flag(examples)+     Buildable:   True+  else+     Buildable:   False++Executable fastout+  Main-Is:        FastOut.hs+  Hs-Source-Dirs: examples .+  Build-Depends:  base >= 3 && <5  +  if flag(examples)+     Buildable:   True+  else+     Buildable:   False++Executable     frecover+  Main-Is:        FRecover.hs+  Hs-Source-Dirs: examples .+  Build-Depends:  base >= 3 && < 5+  Ghc-Options:    -Wall+  if flag(examples)+     Buildable:   True+  else+     Buildable:   False++Executable     frename+  Main-Is:        FRename.hs+  Hs-Source-Dirs: examples .+  Build-Depends:  base >= 3 && < 5, bytestring >= 0.9.1+  Ghc-Options:    -Wall+  if flag(examples)+     Buildable:   True+  else+     Buildable:   False

@@ -0,0 +1,13 @@+{-| FRecover reads a possibly corrupt SFF file, and+  tries to extract as much information as possible from it,+  generating a new SFF file in the the process+-}+module Main where+import Bio.Sequence.SFF+import System.Environment (getArgs)++main :: IO ()+main = mapM_ recoverFile =<< getArgs+  where recoverFile f = writeSFF (f++"_recovered") =<< recoverSFF f+  -- perhaps we should use writeSFF' instead, since it goes back+  -- to update the header with number of reads written?

@@ -0,0 +1,29 @@+{-|+  Rename reads in .SFF files to avoid name clashes.+  Apparently, reads with the same name crashes Newbler, and is+  in any case a bad idea.  This ensures uniqueness by appending a serial number to each read name in a set of files.+-}++module Main where++import Bio.Sequence.SFF+import System.Environment (getArgs)+import qualified Data.ByteString.Char8 as B++main :: IO ()+main = do+  fs <- getArgs+  if null fs then putStrLn "Usage: frename file1.sff [file2.sff ...]"+    else renameSFFs fs+         +renameSFFs :: [FilePath] -> IO ()+renameSFFs = go 0 +  where go _ [] = return ()+        go current (f:fs) = do+          (SFF h rs) <- readSFF f+          writeSFF ("r_"++f) (SFF h $ renameFrom current rs)+          go (current+num_reads h) fs+        renameFrom i rs = zipWith update [i..] rs+          where update j r = let h = read_header r+                                 rn = B.concat [read_name h, B.pack "_", B.pack (show j)]+                             in r { read_header = h { name_length = fromIntegral $ B.length rn, read_name = rn }}

@@ -0,0 +1,19 @@+{-| FastOut reads nucleotide sequences in any supported format,+    and outputs fasta, (sanger) fastq, or illumina fastq.++    Usage: fastout {fasta,fastq,illum} input output+-}++module Main where++import Bio.Sequence+import System.Environment (getArgs)++main = do+  [f,inp,out] <- getArgs+  let out_func = case f of +        "fasta" -> writeFasta+        "fastq" -> writeSangerQ+        "illum" -> writeIllumina+        _ -> error "Usage: fastout {fasta,fastq,illum} input output"+  out_func out =<< readNuc inp

@@ -0,0 +1,21 @@+{-| Flx is a simple tool to trim 454 "Titanium" reads down to "FLX" size, +    i.e. from 800 flows to 400 flows.  This is sometimes useful to investigate+    the relatvie benefit of Titanium's longer read lengths, and will also+    serve as a (rather crude) way to improve average sequence quality.++    We're using 'trimFromTo', which clips a 'ReadBlock' down to specific+    base positions, and we use 'flowToBasePos' to find the base corresponding+    to flow number 400.+-}+  +module Main where++import Bio.Sequence.SFF+import System.Environment (getArgs)+import Data.ByteString as B (take)++-- usage: flx input.sff output.sff+main = do +  [inp,out] <- getArgs+  SFF h rs <- readSFF inp+  writeSFF out (SFF (h { flow_length = 400, flow = B.take 400 (flow h)}) [trimFlows 400 r | r <- rs])