diff --git a/Bio/Alignment/Test.hs b/Bio/Alignment/Test.hs new file mode 100644 --- /dev/null +++ b/Bio/Alignment/Test.hs @@ -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 diff --git a/Bio/Clustering/Test.hs b/Bio/Clustering/Test.hs new file mode 100644 --- /dev/null +++ b/Bio/Clustering/Test.hs @@ -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) diff --git a/Bio/GFF3/Test.hs b/Bio/GFF3/Test.hs new file mode 100644 --- /dev/null +++ b/Bio/GFF3/Test.hs @@ -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 diff --git a/Bio/Location/Test.hs b/Bio/Location/Test.hs new file mode 100644 --- /dev/null +++ b/Bio/Location/Test.hs @@ -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 diff --git a/Bio/Sequence.hs b/Bio/Sequence.hs --- a/Bio/Sequence.hs +++ b/Bio/Sequence.hs @@ -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!" diff --git a/Bio/Sequence/Fasta.hs b/Bio/Sequence/Fasta.hs --- a/Bio/Sequence/Fasta.hs +++ b/Bio/Sequence/Fasta.hs @@ -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) diff --git a/Bio/Sequence/SFF.hs b/Bio/Sequence/SFF.hs --- a/Bio/Sequence/SFF.hs +++ b/Bio/Sequence/SFF.hs @@ -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 ( 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 diff --git a/Bio/Sequence/SFF_filters.hs b/Bio/Sequence/SFF_filters.hs --- a/Bio/Sequence/SFF_filters.hs +++ b/Bio/Sequence/SFF_filters.hs @@ -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 + +-} diff --git a/Bio/Sequence/SeqData.hs b/Bio/Sequence/SeqData.hs --- a/Bio/Sequence/SeqData.hs +++ b/Bio/Sequence/SeqData.hs @@ -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 diff --git a/Bio/Sequence/Test.hs b/Bio/Sequence/Test.hs new file mode 100644 --- /dev/null +++ b/Bio/Sequence/Test.hs @@ -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 diff --git a/Bio/Sequence/TwoBit.hs b/Bio/Sequence/TwoBit.hs --- a/Bio/Sequence/TwoBit.hs +++ b/Bio/Sequence/TwoBit.hs @@ -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 diff --git a/Bio/Util/Test.hs b/Bio/Util/Test.hs new file mode 100644 --- /dev/null +++ b/Bio/Util/Test.hs @@ -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 diff --git a/Bio/Util/TestBase.hs b/Bio/Util/TestBase.hs new file mode 100644 --- /dev/null +++ b/Bio/Util/TestBase.hs @@ -0,0 +1,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) diff --git a/Makefile b/Makefile new file mode 100644 --- /dev/null +++ b/Makefile @@ -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 diff --git a/Test.hs b/Test.hs new file mode 100644 --- /dev/null +++ b/Test.hs @@ -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 + diff --git a/bio.cabal b/bio.cabal --- a/bio.cabal +++ b/bio.cabal @@ -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: . -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 diff --git a/examples/FRecover.hs b/examples/FRecover.hs new file mode 100644 --- /dev/null +++ b/examples/FRecover.hs @@ -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? diff --git a/examples/FRename.hs b/examples/FRename.hs new file mode 100644 --- /dev/null +++ b/examples/FRename.hs @@ -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 }} diff --git a/examples/FastOut.hs b/examples/FastOut.hs new file mode 100644 --- /dev/null +++ b/examples/FastOut.hs @@ -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 diff --git a/examples/Flx.hs b/examples/Flx.hs new file mode 100644 --- /dev/null +++ b/examples/Flx.hs @@ -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])