diff --git a/benchmarks/bench.hs b/benchmarks/bench.hs --- a/benchmarks/bench.hs +++ b/benchmarks/bench.hs @@ -2,22 +2,24 @@ import Criterion.Main import System.Random +import Data.Conduit.Zlib (ungzip, multiple, gzip) import qualified Data.ByteString.Char8 as B import Data.Default.Class -import qualified Data.Conduit.List as CL -import Data.Conduit +import Conduit import Control.Monad.Identity import System.IO.Unsafe import AI.Clustering.Hierarchical +import Data.Either import Bio.Data.Fasta +import Bio.Data.Fastq import Bio.Motif import Bio.Motif.Search import Bio.Motif.Alignment import Bio.Seq dna :: DNA Basic -dna = fromBS $ B.pack $ map f $ take 5000 $ randomRs (0, 3) (mkStdGen 2) +dna = fromRight undefined $ fromBS $ B.pack $ map f $ take 5000 $ randomRs (0, 3) (mkStdGen 2) where f :: Int -> Char f x = case x of @@ -44,10 +46,17 @@ motifs = unsafePerformIO $ readFasta' "data/motifs.fasta" main :: IO () -main = defaultMain - [ bench "motif score" $ nf (scores def pwm) dna - , bgroup "TFBS scanning" [ bench "Naive" $ nf (\x -> runIdentity $ findTFBSSlow def pwm x (0.6 * optimalScore def pwm) $$ CL.consume) dna - , bench "look ahead" $ nf (\x -> runIdentity $ findTFBS def pwm x (0.6 * optimalScore def pwm) $$ CL.consume) dna - ] - , bench "motif merge" $ nf (\x -> fmap show $ flatten $ buildTree x) motifs - ] +main = do + content <- fmap B.concat $ runResourceT $ runConduit $ + sourceFile "tests/data/test.fastq.gz" .| multiple ungzip .| sinkList + let fq = B.concat $ replicate 50 content + defaultMain + [ bench "motif score" $ nf (scores def pwm) dna + --, bgroup "TFBS scanning" [ bench "Naive" $ nf (\x -> runIdentity $ findTFBSSlow def pwm x (0.6 * optimalScore def pwm) $$ CL.consume) dna + -- , bench "look ahead" $ nf (\x -> runIdentity $ findTFBS def pwm x (0.6 * optimalScore def pwm) $$ CL.consume) dna + -- ] + --, bench "motif merge" $ nf (\x -> fmap show $ flatten $ buildTree x) motifs + , bench "Read FASTQ Fast" $ nfAppIO (\x -> runConduit $ yield x .| parseFastqC .| sinkList) fq + , bench "Read FASTQ Slow" $ nfAppIO (\x -> runConduit $ yield x .| parseFastqC' .| sinkList) fq + ] + diff --git a/bioinformatics-toolkit.cabal b/bioinformatics-toolkit.cabal --- a/bioinformatics-toolkit.cabal +++ b/bioinformatics-toolkit.cabal @@ -1,5 +1,5 @@ name: bioinformatics-toolkit -version: 0.9.5 +version: 0.10.0 synopsis: A collection of bioinformatics tools description: A collection of bioinformatics tools license: MIT @@ -10,7 +10,7 @@ category: Bio build-type: Simple extra-source-files: README.md -cabal-version: >=1.18 +cabal-version: 1.18 data-files: tests/data/example.bam tests/data/example.bed.gz @@ -21,7 +21,8 @@ tests/data/example_intersect_peaks.bed.gz tests/data/motifs.fasta tests/data/motifs.meme - tests/data/test.fastq + tests/data/test.fastq.gz + tests/data/test_wrap.fastq.gz tests/data/genes.gtf.gz library @@ -60,9 +61,9 @@ build-depends: base >=4.11 && <5.0 - , aeson + , aeson >= 2.0 , aeson-pretty - , attoparsec + , attoparsec >= 0.13.0.0 , bytestring >= 0.10 , bytestring-lexing >= 0.5 , case-insensitive @@ -73,6 +74,7 @@ , data-ordlist , data-default-class , double-conversion + , deepseq , HsHTSLib >= 1.9.2 , http-conduit >= 2.1.8 , hexpat @@ -108,6 +110,7 @@ , bytestring , data-default-class , conduit + , conduit-extra , mtl test-suite tests diff --git a/src/Bio/Data/Bam.hs b/src/Bio/Data/Bam.hs --- a/src/Bio/Data/Bam.hs +++ b/src/Bio/Data/Bam.hs @@ -17,7 +17,6 @@ import Control.Monad (mzero) import Data.List (foldl') import Bio.Data.Bed -import Bio.Data.Bed.Types import Bio.Data.Fastq import Bio.HTS import Conduit diff --git a/src/Bio/Data/Bed/Types.hs b/src/Bio/Data/Bed/Types.hs --- a/src/Bio/Data/Bed/Types.hs +++ b/src/Bio/Data/Bed/Types.hs @@ -1,5 +1,7 @@ {-# LANGUAGE OverloadedStrings #-} {-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE StrictData #-} +{-# LANGUAGE DeriveGeneric, DeriveAnyClass #-} module Bio.Data.Bed.Types ( BEDLike(..) @@ -32,6 +34,8 @@ import qualified Data.HashMap.Strict as M import qualified Data.IntervalMap.Strict as IM import Data.Maybe (fromJust, fromMaybe) +import GHC.Generics (Generic) +import Control.DeepSeq (NFData) import Bio.Utils.Misc (readDouble, readInt) @@ -87,13 +91,13 @@ -- | BED6 format, as described in http://genome.ucsc.edu/FAQ/FAQformat.html#format1.7 data BED = BED - { _bed_chrom :: !B.ByteString - , _bed_chromStart :: !Int - , _bed_chromEnd :: !Int - , _bed_name :: !(Maybe B.ByteString) - , _bed_score :: !(Maybe Int) - , _bed_strand :: !(Maybe Bool) -- ^ True: "+", False: "-" - } deriving (Eq, Show, Read) + { _bed_chrom :: B.ByteString + , _bed_chromStart :: Int + , _bed_chromEnd :: Int + , _bed_name :: Maybe B.ByteString + , _bed_score :: Maybe Int + , _bed_strand :: Maybe Bool -- ^ True: "+", False: "-" + } deriving (Eq, Show, Read, Generic, NFData) instance Ord BED where compare (BED x1 x2 x3 x4 x5 x6) (BED y1 y2 y3 y4 y5 y6) = @@ -149,10 +153,10 @@ -- | BED3 format data BED3 = BED3 - { _bed3_chrom :: !B.ByteString - , _bed3_chrom_start :: !Int - , _bed3_chrom_end :: !Int - } deriving (Eq, Show, Read) + { _bed3_chrom :: B.ByteString + , _bed3_chrom_start :: Int + , _bed3_chrom_end :: Int + } deriving (Eq, Show, Read, Generic, NFData) instance Ord BED3 where compare (BED3 x1 x2 x3) (BED3 y1 y2 y3) = compare (x1,x2,x3) (y1,y2,y3) @@ -179,11 +183,11 @@ -- | Bedgraph format. data BEDGraph = BEDGraph - { _bdg_chrom :: !B.ByteString - , _bdg_chrom_start :: !Int - , _bdg_chrom_end :: !Int - , _bdg_value :: !Double - } deriving (Eq, Show, Read) + { _bdg_chrom :: B.ByteString + , _bdg_chrom_start :: Int + , _bdg_chrom_end :: Int + , _bdg_value :: Double + } deriving (Eq, Show, Read, Generic, NFData) makeLensesFor [("_bdg_value", "bdgValue")] ''BEDGraph @@ -215,17 +219,17 @@ -- | ENCODE narrowPeak format: https://genome.ucsc.edu/FAQ/FAQformat.html#format12 data NarrowPeak = NarrowPeak - { _npChrom :: !B.ByteString - , _npStart :: !Int - , _npEnd :: !Int - , _npName :: !(Maybe B.ByteString) - , _npScore :: !Int - , _npStrand :: !(Maybe Bool) - , _npSignal :: !Double - , _npPvalue :: !(Maybe Double) - , _npQvalue :: !(Maybe Double) - , _npPeak :: !(Maybe Int) - } deriving (Eq, Show, Read) + { _npChrom :: B.ByteString + , _npStart :: Int + , _npEnd :: Int + , _npName :: Maybe B.ByteString + , _npScore :: Int + , _npStrand :: Maybe Bool + , _npSignal :: Double + , _npPvalue :: Maybe Double + , _npQvalue :: Maybe Double + , _npPeak :: Maybe Int + } deriving (Eq, Show, Read, Generic, NFData) makeLensesFor [ ("_npSignal", "npSignal") , ("_npPvalue", "npPvalue") @@ -276,16 +280,16 @@ -- | ENCODE broadPeak format: https://genome.ucsc.edu/FAQ/FAQformat.html#format13 data BroadPeak = BroadPeak - { _bpChrom :: !B.ByteString - , _bpStart :: !Int - , _bpEnd :: !Int - , _bpName :: !(Maybe B.ByteString) - , _bpScore :: !Int - , _bpStrand :: !(Maybe Bool) - , _bpSignal :: !Double - , _bpPvalue :: !(Maybe Double) - , _bpQvalue :: !(Maybe Double) - } deriving (Eq, Show, Read) + { _bpChrom :: B.ByteString + , _bpStart :: Int + , _bpEnd :: Int + , _bpName :: Maybe B.ByteString + , _bpScore :: Int + , _bpStrand :: Maybe Bool + , _bpSignal :: Double + , _bpPvalue :: Maybe Double + , _bpQvalue :: Maybe Double + } deriving (Eq, Show, Read, Generic, NFData) makeLensesFor [ ("_bpSignal", "bpSignal") , ("_bpPvalue", "bpPvalue") @@ -333,7 +337,7 @@ data BEDExt bed a = BEDExt { _ext_bed :: bed , _ext_data :: a - } deriving (Eq, Show, Read) + } deriving (Eq, Show, Read, Generic, NFData) makeLensesFor [("_ext_bed", "_bed"), ("_ext_data", "_data")] ''BEDExt @@ -346,7 +350,7 @@ strand = _bed . strand instance (Read a, Show a, BEDConvert bed) => BEDConvert (BEDExt bed a) where - asBed chr s e = BEDExt (asBed chr s e) $ error "data not available" + asBed _ _ _ = error "Unable to transform arbitrary record to BEDExt" fromLine l = let (a, b) = B.breakEnd (=='\t') l in BEDExt (fromLine $ B.init a) $ read $ B.unpack b diff --git a/src/Bio/Data/Bed/Utils.hs b/src/Bio/Data/Bed/Utils.hs --- a/src/Bio/Data/Bed/Utils.hs +++ b/src/Bio/Data/Bed/Utils.hs @@ -16,6 +16,7 @@ , queryBaseMap , rpkmBed , rpkmSortedBed + , countTagsBed , countTagsBinBed , countTagsBinBed' , tagCountDistr @@ -176,11 +177,11 @@ V.unsafeFreeze v' let (idx, sortedRegions) = V.unzip v n = G.length idx - rc <- rpkmSortedBed $ Sorted sortedRegions + readCount <- rpkmSortedBed $ Sorted sortedRegions lift $ do result <- GM.new n - G.sequence_ . G.imap (\x i -> GM.unsafeWrite result i (rc U.! x)) $ idx + G.sequence_ . G.imap (\x i -> GM.unsafeWrite result i (readCount U.! x)) $ idx G.unsafeFreeze result {-# INLINE rpkmBed #-} @@ -210,6 +211,27 @@ l = G.length regions {-# INLINE rpkmSortedBed #-} +countTagsBed :: (PrimMonad m, BEDLike b, G.Vector v Int) + => [b] -> ConduitT BED o m (v Int, Int) +countTagsBed regions = do + vec <- lift $ GM.replicate l 0 + n <- foldMC (count vec) (0 :: Int) + vec' <- lift $ G.unsafeFreeze vec + return (vec', n) + where + count v nTags tag = do + let p | tag^.strand == Just True = tag^.chromStart + | tag^.strand == Just False = tag^.chromEnd - 1 + | otherwise = error "Unkown strand" + xs = concat $ IM.elems $ + IM.containing (M.lookupDefault IM.empty (tag^.chrom) intervalMap) p + addOne v xs + return $ succ nTags + intervalMap = bedToTree (++) $ zip regions $ map return [0..] + addOne v' = mapM_ $ \x -> GM.unsafeRead v' x >>= GM.unsafeWrite v' x . (+1) + l = length regions +{-# INLINE countTagsBed #-} + -- | divide each region into consecutive bins, and count tags for each bin and -- return the number of all tags. Note: a tag is considered to be overlapped -- with a region only if the starting position of the tag is in the region. For @@ -226,8 +248,8 @@ v <- GM.replicate num 0 return (v, index) nTags <- foldMC (f vs) 0 - rc <- lift $ mapM (G.unsafeFreeze . fst) $ G.toList vs - return (rc, nTags) + readCount <- lift $ mapM (G.unsafeFreeze . fst) $ G.toList vs + return (readCount, nTags) where f vs n bed = do let pos | bed^.strand == Just True = bed^.chromStart @@ -282,8 +304,8 @@ GM.unsafeRead v i >>= GM.unsafeWrite v i . (+1) sink (nTags+1) vs - _ -> do rc <- lift $ mapM (G.unsafeFreeze . fst) $ G.toList vs - return (rc, nTags) + _ -> do readCount <- lift $ mapM (G.unsafeFreeze . fst) $ G.toList vs + return (readCount, nTags) intervalMap = bedToTree (++) $ zip beds $ map return [0..] {-# INLINE countTagsBinBed' #-} diff --git a/src/Bio/Data/Fastq.hs b/src/Bio/Data/Fastq.hs --- a/src/Bio/Data/Fastq.hs +++ b/src/Bio/Data/Fastq.hs @@ -1,5 +1,6 @@ {-# LANGUAGE LambdaCase #-} {-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE DeriveGeneric, DeriveAnyClass #-} module Bio.Data.Fastq ( Fastq(..) , streamFastqGzip @@ -7,6 +8,7 @@ , sinkFastqGzip , sinkFastq , parseFastqC + , parseFastqC' , fastqToByteString , qualitySummary , trimPolyA @@ -19,6 +21,8 @@ import qualified Data.Attoparsec.ByteString as A import Data.Attoparsec.ByteString.Char8 import Data.Conduit.Attoparsec +import GHC.Generics (Generic) +import Control.DeepSeq -- | A FASTQ file normally uses four lines per sequence. -- @@ -36,7 +40,7 @@ { fastqSeqId :: B.ByteString , fastqSeq :: B.ByteString , fastqSeqQual :: B.ByteString - } deriving (Show, Eq) + } deriving (Show, Eq, Generic, NFData) -- | Read gzipped fastq file. streamFastqGzip :: (PrimMonad m, MonadThrow m, MonadResource m) @@ -55,29 +59,97 @@ => FilePath -> ConduitT Fastq o m () sinkFastqGzip fl = mapC fastqToByteString .| unlinesAsciiC .| gzip .| sinkFileBS fl -parseFastqC :: MonadThrow m => ConduitT B.ByteString Fastq m () -parseFastqC = conduitParser fastqParser .| mapC snd +data FQBuilder = Init FQBuilder + | FQ1 B.ByteString FQBuilder + | FQ2 B.ByteString FQBuilder + | FQ3 B.ByteString FQBuilder + | Complete + +fqBuilder :: FQBuilder -> Fastq +fqBuilder = go ([], [], []) + where + go acc (Init bldr) = go acc bldr + go (f1,f2,f3) (FQ1 x bldr) = go (x:f1, f2, f3) bldr + go (f1,f2,f3) (FQ2 x bldr) = go (f1, x:f2, f3) bldr + go (f1,f2,f3) (FQ3 x bldr) = go (f1, f2, x:f3) bldr + go (f1,f2,f3) Complete = Fastq (B.concat $ reverse f1) + (B.concat $ reverse f2) (B.concat $ reverse f3) +{-# INLINE fqBuilder #-} + +parseFastqC :: Monad m => ConduitT B.ByteString Fastq m () +parseFastqC = await >>= maybe (error "Empty input") ( \x -> do + if B.head x == '@' + then loop Init 'a' $ B.tail x + else error "Record does not start with \'@\'" ) + where + tryRead1 input | B.null input = await >>= maybe (error "Unexpected EOF") return + | otherwise = return input + loop acc st input = case st of + 'a' -> do + (x, rest) <- B.break (=='\n') <$> tryRead1 input + if B.null rest + then loop (acc . FQ1 x) 'a' rest + else loop (acc . FQ1 x) 'b' $ B.tail rest + 'b' -> do + (x, rest) <- B.break (=='\n') <$> tryRead1 input + if B.null rest + then loop (acc . FQ2 x) 'b' rest + else loop (acc . FQ2 x) 'B' $ B.tail rest + 'B' -> do + input' <- tryRead1 input + if B.head input' == '+' + then loop acc 'c' $ B.tail input' + else do + let (x, rest) = B.break (=='\n') input' + if B.null rest + then loop (acc . FQ2 x) 'b' rest + else loop (acc . FQ2 x) 'B' $ B.tail rest + 'c' -> do + (x, rest) <- B.break (=='\n') <$> tryRead1 input + if B.null rest + then loop acc 'c' rest + else loop acc 'd' $ B.tail rest + 'd' -> do + (x, rest) <- B.break (=='\n') <$> tryRead1 input + if B.null rest + then loop (acc . FQ3 x) 'd' rest + else loop (acc . FQ3 x) 'D' $ B.tail rest + 'D' -> if B.null input + then await >>= \case + Nothing -> yield $ fqBuilder $ acc Complete + Just input' -> if B.head input' == '@' + then do + yield $ fqBuilder $ acc Complete + loop Init 'a' $ B.tail input' + else loop acc 'd' input' + else if B.head input == '@' + then do + yield $ fqBuilder $ acc Complete + loop Init 'a' $ B.tail input + else loop acc 'd' input {-# INLINE parseFastqC #-} -fastqParser :: Parser Fastq -fastqParser = do - _ <- skipWhile (/='@') >> char8 '@' - ident <- A.takeTill isEndOfLine - endOfLine - sequ <- BS.filter (not . isEndOfLine) <$> takeTill (=='+') - char8 '+' >> A.skipWhile (not . isEndOfLine) >> endOfLine - score <- BS.filter (not . isEndOfLine) <$> - A.scan 0 (f (B.length sequ)) - skipWhile (/='@') - return $ Fastq ident sequ score +parseFastqC' :: MonadThrow m => ConduitT B.ByteString Fastq m () +parseFastqC' = conduitParser fastqParser .| mapC snd where - f n i x | i >= n = Nothing - | isEndOfLine x = Just i - | otherwise = Just $ i + 1 -{-# INLINE fastqParser #-} + fastqParser = do + _ <- skipWhile (/='@') >> char8 '@' + ident <- A.takeTill isEndOfLine + endOfLine + sequ <- BS.filter (not . isEndOfLine) <$> takeTill (=='+') + char8 '+' >> A.skipWhile (not . isEndOfLine) >> endOfLine + score <- BS.filter (not . isEndOfLine) <$> + A.scan 0 (f (B.length sequ)) + skipWhile (/='@') + return $ Fastq ident sequ score + where + f n i x | i >= n = Nothing + | isEndOfLine x = Just i + | otherwise = Just $ i + 1 +{-# INLINE parseFastqC' #-} fastqToByteString :: Fastq -> B.ByteString -fastqToByteString (Fastq a b c) = "@" <> a <> "\n" <> b <> "\n+\n" <> c +fastqToByteString (Fastq a b c) = B.concat ["@", a, "\n", b, "\n+\n", c] {-# INLINE fastqToByteString #-} -- | Get the mean and variance of quality scores at every position. @@ -87,7 +159,7 @@ meanVarianceC :: Monad m => ConduitT [Double] o m [(Double, Double)] meanVarianceC = peekC >>= \case Nothing -> error "Empty input" - Just x -> fst <$> foldlC f (replicate (length x) (0,0), 0) + Just x -> fst <$> foldlC f (replicate (length x) (0,0), 0 :: Int) where f (acc, n) xs = let acc' = zipWith g acc xs in (acc', n') where diff --git a/src/Bio/GO.hs b/src/Bio/GO.hs --- a/src/Bio/GO.hs +++ b/src/Bio/GO.hs @@ -21,8 +21,6 @@ type GOMap = M.HashMap GOId GO -type TermCount = M.HashMap GOId Int - -- | The top level is 0. getGOLevel :: GOId -> GOMap -> Int getGOLevel gid gm = loop 0 [gid] diff --git a/src/Bio/RealWorld/ENCODE.hs b/src/Bio/RealWorld/ENCODE.hs --- a/src/Bio/RealWorld/ENCODE.hs +++ b/src/Bio/RealWorld/ENCODE.hs @@ -29,12 +29,14 @@ import Data.Aeson import Data.Aeson.Types import Data.Aeson.Encode.Pretty (encodePretty) -import qualified Data.HashMap.Lazy as M +import qualified Data.Aeson.KeyMap as M +import qualified Data.Aeson.Key as M import qualified Data.ByteString.Lazy.Char8 as B import qualified Data.ByteString.Char8 as BS import qualified Data.Sequence as S import qualified Data.Text as T import qualified Data.Vector as V +import Data.Maybe (fromMaybe) import Network.HTTP.Conduit import Data.Default.Class @@ -115,7 +117,7 @@ (|@) :: Value -> T.Text -> Value -(|@) (Object obj) key = M.lookupDefault (error errMsg) key obj +(|@) (Object obj) key = fromMaybe (error errMsg) $ M.lookup (M.fromText key) obj where errMsg = "No such key: " ++ T.unpack key ++ " In: " ++ show obj (|@) _ _ = error "not an object" diff --git a/src/Bio/RealWorld/Ensembl.hs b/src/Bio/RealWorld/Ensembl.hs --- a/src/Bio/RealWorld/Ensembl.hs +++ b/src/Bio/RealWorld/Ensembl.hs @@ -7,7 +7,7 @@ import Data.Aeson import Data.List.Split (chunksOf) import qualified Data.ByteString.Char8 as B -import qualified Data.HashMap.Strict as M +import qualified Data.Aeson.KeyMap as M import Network.HTTP.Conduit import Bio.RealWorld.ID (BioID(..), EnsemblID) diff --git a/src/Bio/Utils/Functions.hs b/src/Bio/Utils/Functions.hs --- a/src/Bio/Utils/Functions.hs +++ b/src/Bio/Utils/Functions.hs @@ -8,7 +8,7 @@ , filterFDR , slideAverage , hyperquick - , kdeWeight + , gaussianKDE , kld , jsd , binarySearch @@ -95,19 +95,19 @@ e = 1e-20 -- | Assign weights to the points according to density estimation. -kdeWeight :: Int -- ^ number of mesh points used in KDE - -> U.Vector Double -> U.Vector Double -kdeWeight n xs = U.map getDensity xs +gaussianKDE :: Int -- ^ number of mesh points used in KDE + -> U.Vector Double -> (Double -> Double) +gaussianKDE n xs = \x -> + let i = binarySearch points x + lo = points U.! (i - 1) + lo_d = den U.! (i - 1) + hi = points U.! i + hi_d = den U.! i + hi_w = (x - lo) / (hi - lo) + lo_w = 1 - hi_w + in lo_w * lo_d + hi_w * hi_d where (points, den) = kde n xs - getDensity x = let i = binarySearch points x - lo = points U.! (i - 1) - lo_d = den U.! (i - 1) - hi = points U.! i - hi_d = den U.! i - hi_w = (x - lo) / (hi - lo) - lo_w = 1 - hi_w - in lo_w * lo_d + hi_w * hi_d -- | compute the Kullback-Leibler divergence between two valid (not check) probability distributions. -- kl(X,Y) = \sum_i P(x_i) log_2(P(x_i)\/P(y_i)). diff --git a/tests/Tests/Fastq.hs b/tests/Tests/Fastq.hs --- a/tests/Tests/Fastq.hs +++ b/tests/Tests/Fastq.hs @@ -3,6 +3,9 @@ module Tests.Fastq (tests) where import Bio.Data.Fastq +import Data.Conduit.Zlib (ungzip, multiple, gzip) +import qualified Data.ByteString.Char8 as B +import qualified Data.ByteString.Lazy.Char8 as BL import Lens.Micro import Conduit import Test.Tasty @@ -14,12 +17,19 @@ tests :: TestTree tests = testGroup "Test: Bio.Data.Fastq" - [ fastqIO + [ testCase "FASTQ IO" fastqIO ] -fastqIO :: TestTree -fastqIO = goldenVsFile "FASTQ Read/Write Test" input output io - where - io = runResourceT $ runConduit $ streamFastq input .| sinkFastq output - input = "tests/data/test.fastq" - output = "out.fastq" +fastqIO :: Assertion +fastqIO = do + a <- fmap B.concat $ runResourceT $ runConduit $ + sourceFile "tests/data/test.fastq.gz" .| multiple ungzip .| sinkList + b <- fmap B.concat $ runResourceT $ runConduit $ + sourceFile "tests/data/test_wrap.fastq.gz" .| multiple ungzip .| sinkList + let fq1 = B.concat $ replicate 1000 a + fq2 = B.concat $ replicate 1000 b + r1 = B.concat $ runIdentity $ runConduit $ yield fq1 .| + parseFastqC .| mapC fastqToByteString .| unlinesAsciiC .| sinkList + r2 = B.concat $ runIdentity $ runConduit $ yield fq2 .| + parseFastqC .| mapC fastqToByteString .| unlinesAsciiC .| sinkList + [True, True] @=? [r1 == fq1, r2 == fq1] diff --git a/tests/data/test.fastq b/tests/data/test.fastq deleted file mode 100644 --- a/tests/data/test.fastq +++ /dev/null @@ -1,200 +0,0 @@ -@M00436:256:000000000-B3PGP:1:1101:12997:1870 2:N:0:TAAGACGA -TTTTTCCCTCCTTTTCTCTTCTCCTCTTTTTTTTTTTTTTTTTTTTTTTTT -+ -111>1111B1111313A33A1B11B01BDAB0AA/AAAAA>E>>>E>EE<< -@M00436:256:000000000-B3PGP:1:1101:20142:2157 2:N:0:TAAGGCGA -CATTAGAGCCTGTTATTAATGATTGTAATGTATTTTCTGTATAATTTTACT -+ -11>>>13B1>>BBGFGBGGGGF3BBAGHFHHHHGG3FBFGHFHFGDFGFHG -@M00436:256:000000000-B3PGP:1:1101:15166:2177 2:N:0:TAAGGCGA -TGCCACGTGATTGCTGAGGCTGTAAGAACCACCCTGGGTCCCCGTGGCATG -+ -A11>>1>>1>FFCGGFGGGGGGHHG3BDFGH0FGCCGGGGHGHGFCCGFHH -@M00436:256:000000000-B3PGP:1:1101:16584:2188 2:N:0:TAAGGCGA -CTAGTAAGGCATTGTAAAAGTTACTGGATTTGGTTTAAAAAAAAAAAAAA -+ ->A11>BF11>CBAA3BGGG11BEGAE1FD1GD1000EGGGHCGGGGG?E@ -@M00436:256:000000000-B3PGP:1:1101:15256:2192 2:N:0:TAAGGCGA -CTTTGGAACAACTATAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA -+ -1AA?11CFBFFABG3FGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG -@M00436:256:000000000-B3PGP:1:1101:18301:2193 2:N:0:TAAGGCGA -CACCTACAGAAGGAGCGCTGTGCCCCCTGGAGCTGACAAGAAAGCTGAGGC -+ -111>AFCC1DFFA1FG?FGAEGHHHGGGGABGGHHDG11000BFGHG1BFG -@M00436:256:000000000-B3PGP:1:1101:13721:2196 2:N:0:TAAGGCGA -TAGCCTGAGCGTTTAATTTATTCAGTACCTGTGTTTGTGTGAATGCGGTGT -+ -111>AFFBBFDAECG3FGDGGGHGGHFHBGHHHBBFFGFGFHBGHGGGEAA -@M00436:256:000000000-B3PGP:1:1101:14584:2200 2:N:0:TAAGGCGA -GCAAAAGGACCATTTAAAGTAAAAAAAAAAAAAAAAAAAAAAAATCTCTCC -+ -1>>1>BCAFFAAGFGG3FGG3FG3GG?ECGGGGGC?EGCEGGE@FH2B2BC -@M00436:256:000000000-B3PGP:1:1101:11335:2200 2:N:0:TAAGGCGA -CTGAATGCCATTAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA -+ -1AAAA33@11113BGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG -@M00436:256:000000000-B3PGP:1:1101:13303:2201 2:N:0:TAAGGCGA -TGCTTATACAATGTGAGTGTACTGAGGTAAGGGGAATATGGGAGCAGAGAC -+ -A3>A?FFF5BCFFDGGGFGGGGHHGHGCF5GHA?EGGFGHGFFEEG3FHCG -@M00436:256:000000000-B3PGP:1:1101:13286:2203 2:N:0:TAAGGCGA -TTCCCGAGACCGGCAGGCCTCAGCCCAAAAGAACAAAGGGAGAAAGAACTG -+ -1>111A>111>10AE00AA101110AGHH/0/00B0B00///A011B0011 -@M00436:256:000000000-B3PGP:1:1101:11054:2203 2:N:0:TAAGGCGA -TTCGTACTCTGAGTTGATACCACTGCTTCCGCGGACAGGCGTGTAGCACTC -+ -1111>>1?DBF1BB31FFE3BA1B1E3BF1A0AEECAF////A/BA22BGC -@M00436:256:000000000-B3PGP:1:1101:15278:2204 2:N:0:TAAGGCGA -ACATAGTTCAGTATGAAAATAAGAGAATAAAATCTGTTATAAGCAAGTGAT -+ -111>1@DFDBFFFGGGGF1BFGF1E11BDHFCGFDHGHHHH3GFH0GAEHH -@M00436:256:000000000-B3PGP:1:1101:14285:2204 2:N:0:TAAGGCGA -AATGTGCGTATTGAGCACATTAAGCACTCTAAGAGCCGAGATAGCTTCCTG -+ -11>1AFBAAA1@D3GGGGGGGGG1GFFH0E1F31E0F?0AECBADDGGHHE -@M00436:256:000000000-B3PGP:1:1101:15307:2205 2:N:0:TAAGGCGA -TCTATAAGAAGGAGAAAGCTCGAGTCATCACTGAGGAAGAGAAGAATTTC -+ -AA3AAFDDBDCFAGECGFFGGGGGHCFHHGHBGHHGHHHHHGHHHHHHBG -@M00436:256:000000000-B3PGP:1:1101:17837:2208 2:N:0:TAAGGCGA -GCAGAGGTGAATGGGATTCAGACGGTACTAAAAAAAAAAAAAAAAAAAAA -+ -111>1A?1>FFFDCFGCEGGGGGCGEHEHFFHFFAEEGCEEGGGGGGGGG -@M00436:256:000000000-B3PGP:1:1101:14152:2208 2:N:0:TAAGGCGA -TCTATGCTCCCTCTGTGTTAGAATTGGCTGAAGTGTTTTTGTGCAGCTGAC -+ -AA>A1BD311>1A113331BFGH3311010BFE313300000A1E1FE1F0 -@M00436:256:000000000-B3PGP:1:1101:17005:2214 2:N:0:TAAGGCGA -GTCTTGAGAGCTTCCTTTTGCTTTCTTCTCTTGCCTCTGAGAACATCAGCC -+ -1AAA>FFCA1A1B1EA3311B313A333313311B11B1F1FGCG1FH110 -@M00436:256:000000000-B3PGP:1:1101:15244:2216 2:N:0:TAAGGCGA -GCCAAGAATGTGTTGTCTAAAATGCCTGTTTAGTTTTCAAGGATGGAACTC -+ -1>>AAAFFD@333B11B3BGGGFAAEA11B3GEAD333FHB1F1EFHH11A -@M00436:256:000000000-B3PGP:1:1101:12952:2217 2:N:0:TAAGGCGA -GAACTGGAGCGTGTGGCCAAGTCTAACCGCTGATTTCCCAGCTGCTGCCT -+ -111>AD1CFBDAEA11FFFGAAGFDFGF0EECE1GFHHFGFHHHFHHGFG -@M00436:256:000000000-B3PGP:1:1101:16915:2218 2:N:0:TAAGGCGA -ATGCACATGTGGTGTAGGAACTTTTTTCTACCATAAGTGACACCAATAAAT -+ -3>3AAFFFFD@BAFFDFGGGGGHHHHGHGHCHGHHHBFBGGFGGGFHGHHH -@M00436:256:000000000-B3PGP:1:1101:11613:2220 2:N:0:TAAGGCGA -TGTGTACTCTGCGTTGATACCACTGCTTCCGCGGACAGGCGTGTAGATCTC -+ -33>AADFFFFFFCGGGGFGGFCHHH4GFHFGDGGGCFFDGGEE1BFFHFH5 -@M00436:256:000000000-B3PGP:1:1101:18615:2225 2:N:0:TAAGGCGA -CCTATCAGGAGTTCACTGACCACCTCGTCAAGACCCACACCAGAGTCTCCG -+ ->AAAAFFFFFB4FFGGFFGGGGHG4GDHGHGBGHBGHGGGGHGFEDGHHHG -@M00436:256:000000000-B3PGP:1:1101:11230:2227 2:N:0:TAAGGCGA -GTACTGTCACTTAACCCCTATTAACATACGGTGTTCAAGCCTTCCAGTATC -+ -?A?AAFFFBDFBFDFFGCGGGGCFHFBDHFG4FE2DFFGHDGBGFGHBGBF -@M00436:256:000000000-B3PGP:1:1101:14569:2228 2:N:0:TAAGGCGA -GTGAATGGAAATGCCCCAACTAAATACTACCGTATGGCCCACCATAATTAC -+ -3>AAAFFFBFFFGFGGFEFFGBGHGFHHHHHGGGGFGHHHGGGDH3GGHGH -@M00436:256:000000000-B3PGP:1:1101:15926:2228 2:N:0:TAAGGCGA -AATCTAGGAAAATAAAACAAGAGAAGCTGAATTAAAAAAAAAAAAAAAAA -+ ->AAA1FF1BDFFGGGGGFGG111B101ABGHDGGHHHHGGGGGGGGGGGG -@M00436:256:000000000-B3PGP:1:1101:15169:2228 2:N:0:TAAGGCGA -GTATCAACGAAAAGTGAGACTCATTGCTGCCCGCCGAAATGAACGGCTGCG -+ -1>A1AFF1>1>AGB11313A11B3331311000A00/A//11B10////// -@M00436:256:000000000-B3PGP:1:1101:15774:2230 2:N:0:TAAGGCGA -CCCAAAGACAAGGAAAGGAAAAAATAAAAACAACAAATACCTAAAAGCTG -+ -111>AFF11>C11111111BE0F0011D110000BGC1F1111BGGD011 -@M00436:256:000000000-B3PGP:1:1101:11431:2230 2:N:0:TAAGGCGA -GAGCTGAACCTGTAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA -+ -111>A1DD1>>1BGGDGGGCEEEEGGGCGGGGGGCEGCEEGCGCGCCCCG -@M00436:256:000000000-B3PGP:1:1101:19030:2230 2:N:0:TAAGGCGA -ATCCCATTACTCTCTTTGTGGAGAAGGAACGAGATAAGGAAGTCAGTGATG -+ -1AAAAFFFFBBD3FFAFF3BEGGHHH00FG?EGGGGFFFHHG1FGHDGGHB -@M00436:256:000000000-B3PGP:1:1101:16869:2231 2:N:0:TAAGGCGA -CCTTCCACCTGGTCCGGAGATCCCTGGACCTTAGGAGGCACTCATTTTTCC -+ -AAAA1>@D111111A1110A01101100A001BA110A/BC0AG12AA10A -@M00436:256:000000000-B3PGP:1:1101:12374:2231 2:N:0:TAAGGCGA -CTTTAGATGGTGGGGACATCATCAATGCCTTGTGCTTCAGCCCCAACCGCT -+ -?AAAAFFFFBBDEEGGEGGGFGHHHG5FGGGHFEDGHHHHGHGGGEGG?EA -@M00436:256:000000000-B3PGP:1:1101:14609:2232 2:N:0:TAAGGCGA -TCCTCCTAAGTAACTGGGACTACATATGTGCATCACCACGTCCAGTTAATT -+ -3>>AAFFFFFFFFGGGFAAG4GG4FD5A66DGHBHH4DC22BFE3GH5EFG -@M00436:256:000000000-B3PGP:1:1101:11898:2233 2:N:0:TAAGGCGA -GTACTGTAGTCACTCAGTATTTGTATATGTTGCTAGAATTTAGATTGTAA -+ -1A1AAFFFFFFFGBGGGBEF3BF3FFGFFGGGHFBGGFFDGGFFFFHFGF -@M00436:256:000000000-B3PGP:1:1101:15215:2233 2:N:0:TAAGGCGA -GTAAAAGAGTATATGTAAACTGAATGAGAACTAAAAAATAGAAATAAAAC -+ -AA1AAFF113@33333BFF1311B3311DF11FGGHH00B11DG1FAGH1 -@M00436:256:000000000-B3PGP:1:1101:18320:2233 2:N:0:TAAGGCGA -GCATTTAATATTGATACACTTTCACCTCTAAAATGGATTTGAAGAAATGC -+ -111>>DF3BBDFGGGGGG1FBGFFCHHHHGF3ADF1GGHFH2A1AB0BFG -@M00436:256:000000000-B3PGP:1:1101:13822:2234 2:N:0:TAAGGCGA -TCCACAAGTACTCTGGCAAGGAAGGAGACAAGCACACCCTGAGCAAGAAGG -+ -BA3ABFA?DFFFFDGGCGGGCGHFE4AGFHBCGFHGFE2EFFFBHHEHBEF -@M00436:256:000000000-B3PGP:1:1101:12165:2235 2:N:0:TAAGGCGA -ACCATGCTGTGTCTGTTATGCGGCAAAGCCGAAGATACTAACAGTATTCTC -+ -AA?A@1B11313A3133B311110AEH0000AE0E1F11F21B11D2222A -@M00436:256:000000000-B3PGP:1:1101:19900:2236 2:N:0:TAAGGCGA -CCATTTTAGTTCAGGAGTGTTTCCTTATTTTCTCTTTTCTGATAGAGTTT -+ ->>3>AFFFFFFFGGGGCGGABGHHHCFHHGGHHHHHHHHHHHHHGHFHHH -@M00436:256:000000000-B3PGP:1:1101:11136:2236 2:N:0:TAAGGCGA -GGAGAGGAGCGAGGCTGAGAGAAGAGAAGAACAATGTACTCATAAAGAAAC -+ ->111111B111110000A0B1AE0F0F101B00B112D2211112111B00 -@M00436:256:000000000-B3PGP:1:1101:15661:2236 2:N:0:TAAGGCGA -GTGTTGCCCAGGCTGGTCTCAAACTCCTGGGCTCAAGCGATCCTCCCACCT -+ ->A>AABFFFFAAGGGFAFFFFGDEGECGF4FDFHFHHHGGEEEHGGG2FEG -@M00436:256:000000000-B3PGP:1:1101:20014:2238 2:N:0:TAAGGCGA -TGTGCATCACAGAGGTGTGCATCACCTTTGACCAGGCTGACCTGACCATC -+ -1111>B3311DDGCC1A1ABGFFHHB1BAGHHHHFEC0FF001FGBHFHH -@M00436:256:000000000-B3PGP:1:1101:13789:2238 2:N:0:TAAGGCGA -ACCCTGGACTCTGGACTTCGCAGGTTCCTGCCTGTCACGCCACCCCCTTCC -+ ->A1A11AFFFBF31FGGCGGGGGGCFFHGBGHGHHHHHGHGGGHGGEEGFH -@M00436:256:000000000-B3PGP:1:1101:11741:2240 2:N:0:TAAGGCGA -CTCCATCCGGCACCCTGACGTTGAGGTGGATGGGTTCTCTGAGCTTCGGTG -+ -1>1A1FDF11DDEECGGFG?FCGHC0B0AFEGCCGCAGBFGGGGGGFGGFE -@M00436:256:000000000-B3PGP:1:1101:19980:2240 2:N:0:TAAGGCGA -GATGTGCGCATTGACACCAGGCTCAACAAAGCTGTCTGGGCCAAAGGAAT -+ -11>1>DBA@AAAGGFFGFGFGGCGHHHHG0FGFHGGHF00/GE/GFFFFH -@M00436:256:000000000-B3PGP:1:1101:20675:2240 2:N:0:TAAGGCGA -GGGTCCCTCTGTACTCTGCGTTGATACCACTGCTTCCGCGGACAGGCGTGT -+ -333AAAFFFFFFGGGGGGCFFGGFGHFHHHHHHHGHHGGGGEFGGGGGEGE -@M00436:256:000000000-B3PGP:1:1101:14624:2241 2:N:0:TAAGGCGA -TACTTATAAGTCTATCATTTAAAGACATGTACTGAAACAAATGTATTTGTT -+ -A1>AAFFFBDDBG3EFGGGGB3FBGHHHFHHHHHFFGFB0BFGEBGGHHGH -@M00436:256:000000000-B3PGP:1:1101:11465:2241 2:N:0:TAAGGCGA -GCCTCACCGGTGACCAGCCTGGACTCGCATGGCAGCACCCAGGGCATGCTG -+ -333AACAFAD>DGGGGGFGGGGHHGHGGEGCHFBEGCHHGGHGG1FGFHHG -@M00436:256:000000000-B3PGP:1:1101:20263:2242 2:N:0:TAAGGCGA -GTGTGAATGGGTTTATACTGATTTTGATATCATGTTCTTTTCATAGTCGT -+ -111>ADDFFC?CCFGGGGGFFGFHHHGHHHCHHHHFHCHHHGHHHHHGFG -@M00436:256:000000000-B3PGP:1:1101:16263:2243 2:N:0:TAAGGCGA -GAACTATACTGACGAAGCCATTGAAACAGATGATCTGACCATTAAGCTGGT -+ -11>AAFFFFFDFGGGG?GGGGGGGFHHHHHGHHBGHGFBFGGFHGFHFF00 diff --git a/tests/data/test.fastq.gz b/tests/data/test.fastq.gz new file mode 100644 Binary files /dev/null and b/tests/data/test.fastq.gz differ diff --git a/tests/data/test_wrap.fastq.gz b/tests/data/test_wrap.fastq.gz new file mode 100644 Binary files /dev/null and b/tests/data/test_wrap.fastq.gz differ