diff --git a/Biobase/Fasta/Export.hs b/Biobase/Fasta/Export.hs deleted file mode 100644 --- a/Biobase/Fasta/Export.hs +++ /dev/null @@ -1,38 +0,0 @@ --- | Fasta export - -module Biobase.Fasta.Export where -import Biobase.Fasta.Types -import qualified Data.ByteString.Lazy.Char8 as B -import Data.List -import GHC.Int - -instance Show Fasta where - show (Fasta _header _sequence) = - (B.unpack _header) ++ "\n" ++ (B.unpack _sequence) ++ "\n" - -prettyPrintFasta :: Int -> Fasta -> String -prettyPrintFasta number (Fasta _header _sequence) = (B.unpack _header) ++ "\n" ++ (B.unpack sequenceLines) ++ "\n" - where sequenceSlices = breakByteString number _sequence - sequenceLines = B.intercalate (B.pack "\n") sequenceSlices - -prettyByteStringFasta :: Int -> Fasta -> B.ByteString -prettyByteStringFasta number (Fasta _header _sequence) = _header `B.append` bslinebreak `B.append` sequenceLines `B.append` bslinebreak - where sequenceSlices = breakByteString number _sequence - sequenceLines = B.intercalate (B.pack "\n") sequenceSlices - bslinebreak = B.pack "\n" - - -breakByteString :: Int -> B.ByteString -> [B.ByteString] -breakByteString number bs - | B.empty == currentLine = [] - | otherwise = currentLine:(breakByteString number rest) - where (currentLine,rest) = B.splitAt (fromIntToInt64 number) bs - -fromIntToInt64 :: Int -> Int64 -fromIntToInt64 = fromIntegral - -writeFastaFile :: FilePath -> [Fasta] -> IO () -writeFastaFile filePath fastas = do - let fastabs = map (prettyByteStringFasta 80) fastas - let outputbs= B.concat fastabs - B.writeFile filePath outputbs diff --git a/Biobase/Fasta/Streaming.hs b/Biobase/Fasta/Streaming.hs --- a/Biobase/Fasta/Streaming.hs +++ b/Biobase/Fasta/Streaming.hs @@ -1,183 +1,196 @@ + -- | Streaming Fasta handling via the @streaming@ library. -- -- The functions in here should be streaming in constant memory. -- +-- A typical, slightly complicated is this: +-- @ +-- forEach :: forall r . Stream (ByteString m) m r -> m (Stream (Of ()) m r) +-- forEach dna = do +-- -- extract the header, but at most 123 characters, dropping the rest +-- hdr SP.:> dta ← extractHeader (Just 123) dna +-- -- create windows @ws@ of a particular type. Include the prefix, the suffix, and make each window 10 characters long +-- let ws = (streamedWindows True True (Just 10) (SequenceIdentifier hdr) PlusStrand dta :: SP.Stream (SP.Of (BioSequenceWindow "DNA" DNA 0)) m r) +-- -- count the number of characters in @dna@, get the return value, print each window +-- count SP.:> r ← SP.mapM_ (liftIO . print) . bswSeqLength $ SP.copy ws +-- liftIO $ print count +-- liftIO $ putStrLn "" +-- -- yield one vacuous @()@ result, return the remainder @r@ from dna. +-- return $ SP.yield () *> return r +-- @ +-- -- TODO Check if this is actually true with some unit tests. -{-# LANGUAGE UnicodeSyntax #-} -{-# LANGUAGE RankNTypes #-} -{-# LANGUAGE LambdaCase #-} -{-# LANGUAGE PolyKinds #-} -{-# LANGUAGE DataKinds #-} - - module Biobase.Fasta.Streaming ( module Biobase.Fasta.Streaming ) where -import Control.Monad -import Control.Monad.Trans.Resource (runResourceT, ResourceT(..), MonadResource) -import Data.ByteString.Streaming as BSS -import Data.ByteString.Streaming.Char8 as S8 -import Data.ByteString.Streaming.Internal (ByteString(..)) -import Data.Semigroup as SG -import Debug.Trace -import GHC.TypeLits -import Prelude as P +import Control.Lens hiding (Index,Empty, mapped) +import Control.Monad +import Control.Monad.Trans.Resource (runResourceT, ResourceT(..), MonadResource) +import Data.Semigroup as SG +import Debug.Trace +import GHC.Generics (Generic) +import GHC.TypeLits +import Prelude as P import qualified Data.ByteString.Char8 as BS -import qualified Data.ByteString.Lazy.Char8 as B import qualified Streaming.Internal as SI -import Streaming as S -import Streaming.Prelude as SP -import qualified Data.List as L -import Biobase.Types.Index.Type -import Biobase.Fasta.Types +import Streaming as S +import Streaming.ByteString as BSS +import Streaming.ByteString.Char8 as S8 +import Streaming.ByteString.Internal as SBI +import Streaming.Prelude as SP +import Data.ByteString.Streaming.Split -newtype HeaderSize = HeaderSize Int - deriving (Eq,Ord,Show) +import Biobase.Types.BioSequence +import Biobase.Types.Index.Type +import Biobase.Types.Location +import Biobase.Types.Position +import Biobase.Types.Strand -newtype OverlapSize = OverlapSize Int - deriving (Eq,Ord,Show) -newtype CurrentSize = CurrentSize Int - deriving (Eq,Ord,Show) -newtype Header (which ∷ k) = Header { getHeader ∷ BS.ByteString } - deriving (Eq,Ord,Show) +-- | -newtype Overlap (which ∷ k) = Overlap { getOverlap ∷ BS.ByteString } - deriving (Eq,Ord,Show) +streamedFasta :: (Monad m) => ByteStream m r -> Stream (Stream (ByteStream m) m) m r +{-# Inlinable streamedFasta #-} +streamedFasta = S.maps collapseData . streamOfStreamedFasta --- | Current Fasta window, together with the start index (0-based). +-- | Here each individual fasta file will be a stream. +-- +-- TODO Once this works, @streamingFasta@ should be @S.concats . streamOfStreamedFasta@ ... -data Current (which ∷ k) = Current { currentFasta ∷ BS.ByteString, currentStart ∷ Index 0 } - deriving (Eq,Ord,Show) +streamOfStreamedFasta + :: forall m r + . ( Monad m ) + => ByteStream m r + -> Stream (Stream (ByteStream m) m) m r + -- ^ +{-# Inlinable streamOfStreamedFasta #-} +streamOfStreamedFasta = go . S8.lines where + go = \case + SI.Return r -> SI.Return r + SI.Effect m -> SI.Effect (fmap go m) + SI.Step fs -> SI.Step (SI.Step (fmap (fmap go . splitFasta) fs)) --- | Fully stream a fasta file, making sure to never exceed a constant amount --- of memory. The @go@ function yields values of type @a@ down the line for --- continued streaming. +-- | Given a 'Stream (ByteString m) m r' which is a 'Stream' of @lines@, split +-- off the first @Fasta@ entry. + +splitFasta :: (Monad m) => Stream (ByteStream m) m r -> Stream (ByteStream m) m (Stream (ByteStream m) m r) +{-# Inlinable splitFasta #-} +splitFasta = loop False where + loop hdr = \case + SI.Return r -> SI.Return (SI.Return r) + SI.Effect m -> SI.Effect (fmap (loop hdr) m) + SI.Step bs -> case bs of + Empty r -> loop hdr r + Chunk cs xs + | BS.null cs -> loop hdr $ SI.Step xs + | h=='>' || h==';' -> if hdr then SI.Return (SI.Step bs) else SI.Step $ fmap (loop True) bs + | otherwise -> SI.Step $ fmap (loop True) bs + where h = BS.head cs + Go m -> SI.Effect $ fmap ((loop hdr) . SI.Step) m + +-- | Given a stream, roughly like @[BS "Header", BS "Data1", BS "Data2", ...]@ +-- create a stream like @[BS "Header", BS "Data"]@. This means that the +-- resulting stream holds exactly two @ByteString@'s. + +collapseData :: (Monad m) => Stream (ByteStream m) m r -> Stream (ByteStream m) m r +{-# Inlinable collapseData #-} +collapseData = loop where + loop = \case + SI.Return r -> SI.Return r + SI.Effect m -> SI.Effect (fmap loop m) + SI.Step bs -> case bs of + Empty r -> loop r + Chunk cs xs + | BS.null cs -> loop $ SI.Step xs + | h=='>' || h==';' -> SI.Step $ fmap (S.yields . S8.concat) bs + | otherwise -> SI.Step $ fmap loop bs + where h = BS.head cs + Go m -> SI.Effect $ fmap (loop . SI.Step) m + + +-- | "Rechunk" a stream of bytestrings. + +reChunkBS :: (Monad m) => Int -> Stream (ByteStream m) m r -> Stream (ByteStream m) m r +{-# Inlinable reChunkBS #-} +reChunkBS n = splitsByteStringAt n . S8.concat + +-- | Assuming a "rechunked" stream of bytestrings, create sequence windows. + +chunksToWindows :: Monad m => SequenceIdentifier w -> Strand -> Stream (ByteStream m) m r -> Stream (Of (Location w FwdPosition (BioSequence ty))) m r +{-# Inlinable chunksToWindows #-} +chunksToWindows seqId s = SP.map go . SP.drop 1 . SP.scan indexed (BS.empty, 0, 0) (\(bs,i,_) -> (bs,i)) . S.mapsM S8.toStrict where + indexed (_,cur,next) bs = (bs,next,next + BS.length bs) + go (bs,i) + = Location + { _locIdentifier = seqId + , _locPosition = FwdPosition s (Index i) + , _locSequence = BioSequence bs + } + + + +-- | Make it possible to take a fasta stream and produce a stream of +-- 'BioSequenceWindow's. This is a convenience function around +-- 'withSuffix . withPrefix . chunksToWindows . reChunks'. -- --- @ --- r4 = toList . streamingFasta (HeaderSize 2) (OverlapSize 1) (CurrentSize 2) go . S8.fromStrict $ BS.pack t0 --- where go (Header h) (Overlap o) (Current c) = yield (h,o,c) --- @ +-- In case of a @Nothing@ window size, a single huge @Fasta@ entry is produced +-- (and materialized!). +-- +-- TODO In case of @Nothing@ window size, we use the 'collapseData' function +-- which has one check too many, and will be slightly slower. However, the +-- check should be once per @ByteString@. -streamingFasta - ∷ forall m w r a - . ( Monad m ) - ⇒ HeaderSize - -- ^ Maximal length of the header. Ok to set to @20 000@, only guards against - -- an extremely long header line. - → OverlapSize - -- ^ How much of the current size to carry over to the next step. Even if set - -- larger than current size, it will only be at most current size. (But see - -- todo at 'overlappedFasta') - → CurrentSize - -- ^ The size of each window to be processed. - → (Header w → Overlap w → Current w → Stream (Of a) m ()) - -- ^ The processing function. Takes in the header, any overlap from the - -- previous window, the current window and produces a stream of @a@s. - → ByteString m r - -- ^ A streaming bytestring of Fasta files. - → Stream (Of a) m r - -- ^ The outgoing stream of @a@s being processed. -{-# Inlinable streamingFasta #-} -streamingFasta (HeaderSize hSz) (OverlapSize oSz) (CurrentSize cSz) f = go (FindHeader [] 0) where - -- Find the next FASTA header - go (FindHeader hdr cnt) = \case - -- No more data to be had. If There is some part of a header, we will run - -- the handling function @f@ with empty input. @f@ can decide on how to - -- handle empty FASTA entries. - Empty retVal → do - -- handle case of last empty fasta - unless (P.null hdr) $ do - let thisHeader = BS.take hSz $ BS.concat $ P.reverse hdr - f (Header thisHeader) (Overlap BS.empty) (Current BS.empty 0) - SI.Return retVal - -- Effects are wrapped up into a 'Stream' effect. - Go m → SI.Effect $ liftM (go (FindHeader hdr cnt)) m - -- We have a chunk of bytestring @rawBS@ with more data in the bytestream - -- @bs@. We work on @b@, not the @rawBS@. In case we have no header parts - -- yet, all characters preceeding a fasta header symbol ('>' or ';') are - -- dropped. - Chunk rawBS bytestream - -- No newline in the @b@, hence we add the bytestring to the partial - -- header, and continue scanning. Note that we add only if we are below - -- the maximal header size @hSz@ to prevent malicious fasta files from - -- blowing up memory usage. - | Nothing ← mk → if cnt > hSz - then go (FindHeader hdr cnt) bytestream - else go (FindHeader (b:hdr) (BS.length b + cnt)) bytestream - -- We have found a newline at @k@. Prepare the full header (up to @hSz@ - -- size) and hand over to @HasHeader@ which processes actual fasta - -- payload. - | Just k ← mk → let thisHeader = BS.take hSz $ BS.concat $ P.reverse $ BS.take k b:hdr - in go (HasHeader thisHeader BS.empty [] 0 0) - (Chunk (BS.drop (k+1) b) bytestream) - where b = if P.null hdr then BS.dropWhile (\c → c/='>' && c/=';') rawBS else rawBS - mk = BS.elemIndex '\n' b - -- We actually do have a valid header now and process fasta in parts. - go hasHeader@(HasHeader hdr overlap cs cnt entries) = \case - -- No more data, process final input and return. - Empty retVal → do - when (cnt>0 || entries==0) $ f (Header hdr) (Overlap BS.empty) (Current (BS.concat $ reverse cs) 0) - SI.Return retVal - -- Effects to be dealt with. - Go m → SI.Effect $ liftM (go hasHeader) m - -- We have incoming data ... - Chunk b bytestream → case newFastaIndex b of - -- there is no new fasta starting, meaning that we need to process @b@ as - -- payload. We split at the maximal size we are allowed according to - -- @cSz@. If we have hit the limit, we run @f@ on this part of the data - -- and include the overlap as prefix. Otherwise we continue gathering. - -- Any newlines are removed from the data. - Nothing → let (this,next) = BS.splitAt (cSz-cnt) $ BS.filter (/= '\n') b - in if BS.length this + cnt >= cSz - then do let thisFasta = BS.concat $ reverse $ this:cs - f (Header hdr) (Overlap overlap) (Current thisFasta 0) - go (HasHeader hdr (BS.drop (BS.length thisFasta - oSz) thisFasta) [] 0 (entries+1)) - (if BS.null next then bytestream else Chunk next bytestream) - else go (HasHeader hdr overlap (this:cs) (BS.length this + cnt) entries) - (if BS.null next then bytestream else Chunk next bytestream) - -- We have a new fasta symbol in @b@. We split at the symbol and re-run - -- the first part (which will end up being the @Nothing@ case) and put - -- into @Chunk next bytestream@ the beginning of the next fasta entry. - -- This part will then be handled by the @otherwise@ case here. - Just new - | new > 0 → let (this,next) = BS.splitAt new b - in go (HasHeader hdr overlap cs cnt entries) $ Chunk this (Chunk next bytestream) - | otherwise → do let thisFasta = BS.concat $ reverse cs - -- we only emit on empty @thisFasta@, if there is - -- data, or it is the only (then empty) entry. - when (cnt>0 || entries==0) $ f (Header hdr) (Overlap overlap) (Current thisFasta 0) - go (FindHeader [] 0) $ Chunk b bytestream - -- Returns the first index (if any) of a new fasta entry symbol. - newFastaIndex b = getMin <$> (Min <$> BS.elemIndex '>' b) SG.<> (Min <$> BS.elemIndex ';' b) +streamedWindows + :: (Monad m) + => Maybe Int + -> Maybe Int + -> Maybe Int + -- ^ desired size or a single huge @Fasta@ entry. + -> SequenceIdentifier w + -> Strand + -> (Stream (ByteStream m) m) r +-- -> Stream (Of (BioSequenceWindow w ty FwdLocation)) m r + -> Stream (Of (PIS w FwdPosition (BioSequence ty))) m r +{-# Inlinable streamedWindows #-} +streamedWindows withPrefix withSuffix winSz seqId strnd + = (maybe id attachSuffixes withSuffix) + . (maybe id attachPrefixes withPrefix) + . SP.map pis + . chunksToWindows seqId strnd + . (case winSz of { Nothing -> collapseData; Just sz -> reChunkBS sz }) --- | Control structure for 'streamingFasta'. +-- | Get the full length of a stream of 'BioSequenceWindow's, counted in +-- characters in each 'bswSequence'. +-- +-- To use, start with @bswSeqLength $ SP.copy xs@. Then consume this stream +-- normally. It still provides a 'Stream' of 'BioSequenceWindows's. However, +-- the return type is now not just @r@, but it provides @Int SP.:> r@, where +-- the @Int@ provides the total length of characters within this @Fasta@ entry. +-- +-- This value may then be used to fully update negative strand information. -data FindHeader - = FindHeader - { headerParts ∷ [BS.ByteString] - -- ^ the collected header parts (in reverse order) - , headerLength ∷ !Int - -- ^ accumulated header length - } - | HasHeader - { header ∷ !BS.ByteString - -- ^ the (size-truncated) header for this fasta file - , dataOverlap ∷ !BS.ByteString - -- ^ overlap (if any) from earlier parts of the fasta file - , dataParts ∷ [BS.ByteString] - -- ^ collection of dataParts, in reverse order! - , dataLength ∷ !Int - -- ^ total length of data parts, simplifies checking if enough data was collected - , entries ∷ !Int - -- ^ count how many entries we have seen - } +streamLocationLength :: (Monad m, ModifyLocation posTy seqTy) => Stream (Of (Location i posTy seqTy)) m r -> m (Of Int r) +{-# Inlinable streamLocationLength #-} +streamLocationLength = SP.fold (\x w -> x + locLength w) 0 id +-- | As a first function, the header should be extracted from a @Fasta@ stream. Since headers may be +-- malformed / malicious, we make it possible to +extractHeader + :: (Monad m) + => Maybe Int + -> Stream (ByteStream m) m r + -> m (Of BS.ByteString (Stream (ByteStream m) m r)) +{-# Inlinable extractHeader #-} +extractHeader hdrSz = + let go = case hdrSz of { Nothing -> id; Just sz -> S8.drained . S8.splitAt (fromIntegral sz) } + in S8.toStrict . go . S8.concat . S.splitsAt 1 + + {- t0 = P.unlines [ ">Aaaa" @@ -189,40 +202,30 @@ ] -r2 = splitFastaLines $ S8.lines $ S8.fromStrict $ BS.pack t0 - -r3 = streamFastaLines $ S8.lines $ S8.fromStrict $ BS.pack t0 - --- r3' ∷ Stream (Stream (Of BS.ByteString) Identity) Identity () -r3' = toList . mapped toList $ maps (mapped toStrict) r3 - -r4 = toList . streamingFasta (HeaderSize 2) (OverlapSize 1) (CurrentSize 2) go . S8.fromStrict $ BS.pack t0 - where go (Header h) (Overlap o) (Current c) = yield (h,o,c) +r4 = toList . streamingFasta (HeaderSize 2) (OverlapSize 1) (CurrentSize 2) . S8.fromStrict $ BS.pack t0 -} ---eachFasta :: forall (m0 :: * -> *). Header Int -> Overlap Int -> Current Int -> Stream (Of (BS.ByteString, BS.ByteString, BS.ByteString)) (ResourceT IO) () -eachFasta (Header h) (Overlap o) (Current c p) = SP.yield (h,o,c) - ---readFastaFile ∷ FilePath → IO () -- [(BS.ByteString,BS.ByteString,BS.ByteString)] ---readFastaFile f = do --- let s = 1000000000000 --- r ← runResourceT --- $ SP.mapM_ (liftIO . P.print) --- $ streamingFasta (HeaderSize s) (OverlapSize 0) (CurrentSize s) eachFasta --- $ S8.readFile f --- return r +{- +--eachFasta (Header h) (Overlap o) (Current c p) = SP.yield (h,o,c) +eachFasta (Header h) (Overlap o) (Current c p) = SP.yield (BS.length h, BS.length o, BS.length c) -parseFastaFile ∷ FilePath → IO [Fasta] -parseFastaFile f = do +--readFastaFile :: FilePath -> IO [(BS.ByteString,BS.ByteString,BS.ByteString)] +readFastaFile f = do let s = 1000000000000 r ← runResourceT - $ toList_ + $ SP.mapM_ (liftIO . P.print) $ streamingFasta (HeaderSize s) (OverlapSize 0) (CurrentSize s) eachFasta $ S8.readFile f - let fastas = L.map (\(a,_,c) -> Fasta (B.fromStrict a) (B.fromStrict c)) r - return fastas + return r +-} -parseFasta ∷ B.ByteString → [Fasta] -parseFasta input = L.map (\(a,_,c) -> Fasta (B.fromStrict a) (B.fromStrict c)) (L.head r) - where s = 1000000000000 - r = toList_ $ streamingFasta (HeaderSize s) (OverlapSize 0) (CurrentSize s) eachFasta $ BSS.fromLazy input +{- +readFastaFile f = do + let s = 1000000000000 + r ← runResourceT + $ SP.mapM_ (liftIO . P.print) + $ SP.mapped S8.toStrict + $ S8.split '>' + $ S8.readFile f + return r +-} diff --git a/Biobase/Fasta/Strict.hs b/Biobase/Fasta/Strict.hs new file mode 100644 --- /dev/null +++ b/Biobase/Fasta/Strict.hs @@ -0,0 +1,107 @@ + +-- | A convenience module for *small* @Fasta@ entries, that are completely in +-- memory and *not* to be streamed. +-- +-- The @Data.ByteString.Strict.Lens@ module is very helpful for further +-- handling of 'Fasta' entries. +-- +-- For convenience, the 'convertString' function from @string-conversions@ is +-- supplied. + +module Biobase.Fasta.Strict + ( module Biobase.Fasta.Strict + , convertString + ) where + +import Control.Lens +import Data.Bifunctor (first) +import Data.ByteString (ByteString) +import Data.String.Conversions +import Data.Void +import GHC.Generics (Generic) +import qualified Data.ByteString.Builder as BB +import qualified Data.ByteString.Char8 as BS +import qualified Data.ByteString.Lazy as BSL +import qualified Data.ByteString.Streaming as BSS +import qualified Streaming.Prelude as SP + +import Biobase.Fasta.Streaming as FS +import Biobase.Types.BioSequence + + + +-- | A *strict* @Fasta@ entry. + +data Fasta which ty = Fasta + { _header ∷ !(SequenceIdentifier which) + , _fasta ∷ !(BioSequence ty) + } + deriving (Eq,Ord,Read,Show,Generic) +makeLenses ''Fasta + +-- | If you don't want to deal with the phantom types. + +type FastaUntyped = Fasta Void Void + +-- | Render a 'Fasta' entry to a 'ByteString'. Will end with a final @\n@ in +-- any case. + +fastaToByteString ∷ Int → Fasta which ty → ByteString +{-# Inlinable fastaToByteString #-} +fastaToByteString k' Fasta{..} = BS.cons '>' (_header^._Wrapped) <> "\n" <> go (_fasta^._Wrapped) + where go (BS.splitAt k → (hd,tl)) + | BS.null hd = mempty + | otherwise = hd <> "\n" <> go tl + k = max 1 k' + +-- | Render a 'Fasta' entry to a 'Builder'. Will end with a final @\n@ in +-- any case. + +fastaToBuilder ∷ Int → Fasta which ty → BB.Builder +{-# Inlinable fastaToBuilder #-} +fastaToBuilder k' Fasta{..} = BB.char8 '>' <> (BB.byteString $ _header^._Wrapped) <> BB.char8 '\n' <> go (_fasta^._Wrapped) + where go (BS.splitAt k → (hd,tl)) + | BS.null hd = mempty + | otherwise = BB.byteString hd <> BB.char8 '\n' <> go tl + k = max 1 k' + +-- | Try to parse a 'ByteString' as a 'Fasta', failing with 'Left', succees +-- with 'Right'. + +byteStringToFasta ∷ ByteString → Either String (Fasta which ty) +{-# Inlinable byteStringToFasta #-} +byteStringToFasta (BS.lines → ls) + | null ls = Left "empty bytestring" + | Just (z, hdr) ← BS.uncons h, z `BS.elem` ">;" = Right $ Fasta { _header = SequenceIdentifier hdr, _fasta = BioSequence $ BS.concat ts } + | otherwise = Left "no '>'/';' first character" + where h:ts = ls + +-- | Try to parse a 'ByteString' as multiple 'Fasta' entries. Even though this +-- is using the underlying streaming interface, this is not streaming. + +{- +byteStringToMultiFasta + ∷ BSL.ByteString → [Fasta which ty] +{-# Inlinable byteStringToMultiFasta #-} +byteStringToMultiFasta bsl = map (view windowedFasta) $ runIdentity bss + where bss = SP.toList_ . streamingFasta (HeaderSize maxBound) (OverlapSize 0) (CurrentSize maxBound) $ BSS.fromLazy bsl +-} + +-- | A lens that goes from a 'BioSequenceWindow' to a 'Fasta'. + +{- +windowedFasta ∷ Lens' (BioSequenceWindow w ty k) (Fasta w ty) +{-# Inline windowedFasta #-} +windowedFasta = lens lr rl + where lr bsw = Fasta { _header = bsw^.bswIdentifier, _fasta = bsw^.bswSequence } + rl bsw f = set bswSequence (f^.fasta) $ set bswIdentifier (f^.header) bsw +-} + +-- | A prism from a 'ByteString' to a 'Fasta'. Note that this will only be an +-- identity if the underlying fasta file is rendered with @k@ characters per +-- line. + +rawFasta ∷ Int → Prism' ByteString (Fasta which ty) +{-# Inline rawFasta #-} +rawFasta k = prism (fastaToByteString k) $ \bs → first (const bs) $ byteStringToFasta bs + diff --git a/Biobase/Fasta/Types.hs b/Biobase/Fasta/Types.hs deleted file mode 100644 --- a/Biobase/Fasta/Types.hs +++ /dev/null @@ -1,56 +0,0 @@ -{-# Language DeriveGeneric #-} - -module Biobase.Fasta.Types where - -import Control.DeepSeq -import Control.Lens -import Data.ByteString.Char8 (ByteString) -import qualified Data.ByteString.Lazy.Char8 as B -import Data.Data -import GHC.Generics -import Biobase.Types.NucleotideSequence -import Biobase.Types.AminoAcidSequence - --- | - -data Fasta = Fasta { fastaHeader :: B.ByteString, fastaSequence :: B.ByteString } - deriving (Eq) - -newtype RawFastaEntry = RawFastaEntry { _rawFastaEntry :: ByteString } - deriving (Show,Eq,Ord,Typeable) - ---makeLenses ''RawFastaEntry - --- | 'StreamEvent's are chunked pieces of data, where the raw data is --- a strict @ByteString@. Each element also retains information on the --- first and last line and column (via 'streamLines') that are part of this --- chunk. - -data StreamEvent - -- | A Header event, multiple header events signal that the header name - -- was longer than the chunk size. - = StreamHeader { streamHeader :: !ByteString, streamLines :: !LineInfo } - -- | A data event. We keep a pointer to the previous chunk (which is - -- useful for some algorithms). The chunk is free of newlines! - | StreamFasta { streamFasta :: !ByteString, prevStreamFasta :: !ByteString, streamLines :: !LineInfo, streamHeader :: !ByteString } - deriving (Show,Eq,Ord,Typeable,Generic) - -instance NFData StreamEvent - - - --- | Complete information on line and column start and end for a chunk. --- --- TODO This is a 1-based format? Lets use the BiobaseTypes facilities! - -data LineInfo = LineInfo - { firstLine :: !Int -- ^ first line for this chunk @(lines in complete file!)@ - , firstCol :: !Int -- ^ first column in first line for this chunk - , lastLine :: !Int -- ^ last line for this chunk @(lines in complete file!)@ - , lastCol :: !Int -- ^ last column in last line for this chunk - , firstIndex :: !Int -- ^ first index in this fasta block. Counts just the number of symbols in the @Fasta@ payload. - } - deriving (Show,Eq,Ord,Typeable,Generic) - -instance NFData LineInfo - diff --git a/BiobaseFasta.cabal b/BiobaseFasta.cabal --- a/BiobaseFasta.cabal +++ b/BiobaseFasta.cabal @@ -1,17 +1,17 @@ +cabal-version: 2.2 name: BiobaseFasta -version: 0.2.0.0 -author: Christian Hoener zu Siederdissen, Florian Eggenhofer +version: 0.4.0.1 +author: Christian Hoener zu Siederdissen maintainer: choener@bioinf.uni-leipzig.de homepage: https://github.com/choener/BiobaseFasta bug-reports: https://github.com/choener/BiobaseFasta/issues -copyright: Christian Hoener zu Siederdissen, 2011-2018 +copyright: Christian Hoener zu Siederdissen, 2011-2021 category: Bioinformatics -license: GPL-3 +license: BSD-3-Clause license-file: LICENSE build-type: Simple stability: experimental -cabal-version: >= 1.10.0 -tested-with: GHC == 8.4.3 +tested-with: GHC == 8.8.4, GHC == 8.10.4, GHC == 9.0 synopsis: streaming FASTA parser description: Stream-based handling of FASTA files. The user selects a window @@ -19,9 +19,6 @@ previous (past) window is available, in case some data sits on the boundary between windows. . - FastaTool is a simple tool providing information on FASTA - files, and allowing to extract sequences and subsequences. - . Greg Schwartz' package is a lot more complete. This one is mostly tailored to my usage requirements (and may at some point use his library). @@ -37,20 +34,17 @@ -library +common deps build-depends: base >= 4.7 && < 5.0 , bytestring + , lens >= 4.0 , resourcet >= 1.0 , streaming >= 0.1 - , streaming-bytestring >= 0.1 - , lens - , deepseq + , streaming-bytestring >= 0.2 + , string-conversions >= 0.4 -- - , BiobaseTypes == 0.1.4.* - exposed-modules: - Biobase.Fasta.Streaming, - Biobase.Fasta.Types, - Biobase.Fasta.Export + , BiobaseTypes == 0.2.1.* + , DPutils == 0.1.1.* default-language: Haskell2010 default-extensions: BangPatterns @@ -58,24 +52,56 @@ , DeriveDataTypeable , DeriveGeneric , FlexibleContexts + , FlexibleInstances , GADTs + , GeneralizedNewtypeDeriving , KindSignatures , LambdaCase + , MultiParamTypeClasses , MultiWayIf , NoMonomorphismRestriction + , OverloadedStrings , PolyKinds , RankNTypes , RecordWildCards , ScopedTypeVariables , TemplateHaskell + , TypeApplications + , TypeFamilies , UnicodeSyntax , ViewPatterns ghc-options: -O2 +library + import: + deps + exposed-modules: + Biobase.Fasta.Streaming + Biobase.Fasta.Strict + +-- | A simple tool for fasta files, showing some features + +executable fastaextract + import: deps + build-depends: base + , optparse-applicative >= 0.14 + -- + , BiobaseFasta + hs-source-dirs: + src + main-is: + fastaextract.hs + ghc-options: + -rtsopts + + + test-suite properties + import: + deps type: exitcode-stdio-1.0 main-is: @@ -84,19 +110,8 @@ -threaded -rtsopts -with-rtsopts=-N hs-source-dirs: tests - default-language: - Haskell2010 - default-extensions: TemplateHaskell - , UnicodeSyntax - , OverloadedStrings - - build-depends: base - , QuickCheck - , bytestring + build-depends: QuickCheck , filepath - , resourcet >= 1.0 - , streaming >= 0.1 - , streaming-bytestring >= 0.1 , tasty >= 0.11 , tasty-hunit >= 0.9 , tasty-golden >= 2.3 diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -1,4 +1,5 @@ -[![Build Status](https://travis-ci.org/choener/BiobaseFasta.svg?branch=master)](https://travis-ci.org/choener/BiobaseFasta) +![github action: master](https://github.com/choener/BiobaseFasta/actions/workflows/ci.yml/badge.svg?branch=master) +![github action: hackage](https://github.com/choener/BiobaseFasta/actions/workflows/hackage.yml/badge.svg) # BiobaseFasta diff --git a/changelog.md b/changelog.md --- a/changelog.md +++ b/changelog.md @@ -1,3 +1,9 @@ +0.4.0.1 +------- + +- streaming-bytestring >= 0.2 +- cleanup of older code + 0.2.0.0 ------- diff --git a/src/fastaextract.hs b/src/fastaextract.hs new file mode 100644 --- /dev/null +++ b/src/fastaextract.hs @@ -0,0 +1,60 @@ + +module Main where + +import Control.Monad.IO.Class +import Data.ByteString.Char8 as BS +import Data.ByteString.Streaming.Char8 as BSS +import Options.Applicative +import Streaming as S +import Streaming.Prelude as SP +import System.IO (stdin) + +import Biobase.Fasta.Streaming as FS + + + +data Options + -- Extract all sequences that have "header" as "infix". + = Extract + { header ∷ String + , from ∷ Int + , to ∷ Int + } + +options ∷ Parser Options +options + = Extract + <$> strOption (long "header" <> short 'h' <> help "header infix to grep") + <*> option auto (long "from" <> short 'f' <> help "first nucleotide in sequence") + <*> option auto (long "to" <> short 't' <> help "last nucleotide in sequence") + +-- | Extract a fasta piece from a larger fasta + +extract + ∷ ( Monad m ) + ⇒ String → Int → Int + → Stream (BSS.ByteString m) m r + → BSS.ByteString m r +{-# Inlinable extract #-} +extract ifx' f' t' s = BSS.mwrap $ do + let ifx = BS.pack ifx' + let f = fromIntegral $ min f' t' + let t = fromIntegral $ max f' t' + hdr :> dta ← extractHeader Nothing s + if ifx `BS.isInfixOf` hdr + -- we actually have a stream to return + then return $ do + BSS.fromStrict $ hdr `BS.snoc` '\n' + BSS.drained . BSS.splitAt (t-f+1) . BSS.drop (f-1) $ BSS.concat dta + -- just drain this stream + else mapsM_ BSS.effects dta >>= return . return + +main ∷ IO () +main = do + p ← execParser (info options fullDesc) + case p of + Extract hdr f t + → BSS.stdout . BSS.unlines . BSS.denull + . maps (extract hdr f t) + . FS.streamedFasta $ BSS.stdin + diff --git a/tests/properties.hs b/tests/properties.hs --- a/tests/properties.hs +++ b/tests/properties.hs @@ -1,63 +1,70 @@ module Main where -import Prelude as P -import Data.Functor.Of -import qualified Data.ByteString.Char8 as BS +import Control.Lens (view) +import Control.Monad.Trans.Resource (runResourceT, ResourceT(..), MonadResource) import Data.ByteString.Streaming as BSS import Data.ByteString.Streaming.Char8 as S8 import Data.ByteString.Streaming.Internal (ByteString(..)) +import Data.Functor.Of +import Data.Void +import Prelude as P +import qualified Data.ByteString.Char8 as BS +import qualified Data.ByteString.Lazy.Char8 as BL8 +import qualified Data.Text as T +import qualified Data.Text.Encoding as T +import qualified Test.Tasty.Golden as Golden import Streaming as S import Streaming.Prelude as SP +import System.FilePath (takeBaseName, replaceExtension) import Test.Tasty import Test.Tasty.HUnit import Test.Tasty.QuickCheck as QC ---import Test.Tasty.Silver as S ---import Test.Tasty.Silver.Interactive as SI import Test.Tasty.TH -import Control.Monad.Trans.Resource (runResourceT, ResourceT(..), MonadResource) -import qualified Data.Text as T -import qualified Data.Text.Encoding as T -import qualified Test.Tasty.Golden as Golden -import System.FilePath (takeBaseName, replaceExtension) -import qualified Data.ByteString.Lazy.Char8 as BL8 +import Biobase.Types.BioSequence +import Biobase.Types.Location +import Biobase.Types.Strand + import Biobase.Fasta.Streaming +import Biobase.Fasta.Strict --- * golden tests - -readFastaFile ∷ FilePath → IO [(BS.ByteString,BS.ByteString,BS.ByteString)] -readFastaFile f = do - let s = 1000000 - xs :> r ← runResourceT - $ SP.toList - $ streamingFasta (HeaderSize s) (OverlapSize 0) (CurrentSize s) eachFasta - $ S8.readFile f - return xs - -lazyFastaOut = BL8.concat . P.map go - where go (h,o,c) = BL8.concat - [ BL8.fromStrict h - , BL8.pack "\n" - , BL8.fromStrict o - , BL8.pack "\n" - , BL8.fromStrict c - , BL8.pack "\n" - ] - -goldenTests ∷ IO TestTree -goldenTests = do - fastaFiles ← Golden.findByExtension [".fa"] "./tests/" - return $ testGroup "readFastaFile golden tests" - [ Golden.goldenVsString - (takeBaseName fastaFile) -- test name - goldenFile -- golden file path - (lazyFastaOut <$> readFastaFile fastaFile) -- action whose result is tested - | fastaFile <- fastaFiles - , let goldenFile = replaceExtension fastaFile ".fa-golden" - ] +-- -- * golden tests +-- +-- eachFasta (HeaderSize h) (OverlapSize o) (CurrentSize c p) = SP.yield (h,o,c) +-- +-- readFastaFile ∷ FilePath → IO [(BS.ByteString,BS.ByteString,BS.ByteString)] +-- readFastaFile f = do +-- let s = 1000000 +-- xs :> r ← runResourceT +-- $ SP.toList +-- $ streamingFasta (HeaderSize s) (OverlapSize 0) (CurrentSize s) eachFasta +-- $ S8.readFile f +-- return xs +-- +-- lazyFastaOut = BL8.concat . P.map go +-- where go (h,o,c) = BL8.concat +-- [ BL8.fromStrict h +-- , BL8.pack "\n" +-- , BL8.fromStrict o +-- , BL8.pack "\n" +-- , BL8.fromStrict c +-- , BL8.pack "\n" +-- ] +-- +-- goldenTests ∷ IO TestTree +-- goldenTests = do +-- fastaFiles ← Golden.findByExtension [".fa"] "./tests/" +-- return $ testGroup "readFastaFile golden tests" +-- [ Golden.goldenVsString +-- (takeBaseName fastaFile) -- test name +-- goldenFile -- golden file path +-- (lazyFastaOut <$> readFastaFile fastaFile) -- action whose result is tested +-- | fastaFile <- fastaFiles +-- , let goldenFile = replaceExtension fastaFile ".fa-golden" +-- ] -- * unit tests @@ -70,28 +77,35 @@ , "890" ] -smallTest ∷ Int → Int → Int → Of [(BS.ByteString,BS.ByteString,BS.ByteString)] () -smallTest h o c = runIdentity - . toList - . streamingFasta (HeaderSize h) (OverlapSize o) (CurrentSize c) go - . S8.fromStrict - $ BS.pack smallInlineFasta - where go (Header h) (Overlap o) (Current c _) = yield (h,o,c) - -smallTest333 = testCase "3/3/3" $ do - let res :> r = smallTest 3 3 3 - assertEqual "return is null" () r - assertEqual "length is 4" 4 (P.length res) - assertEqual "!!0" (">Aa","","123") (res!!0) - assertEqual "!!1" (">Bb","","456") (res!!1) - assertEqual "!!2" (">Bb","456","7") (res!!2) - assertEqual "!!3" (">Cc","","890") (res!!3) +--smallTest ∷ Int → Int → Int → Of [BioSequenceWindow Void Void PartialLocation] () +--smallTest h o c = runIdentity +-- . toList +---- . SP.map (view windowedFasta) +-- . streamingFasta (HeaderSize h) (OverlapSize o) (CurrentSize c) +-- . S8.fromStrict +-- $ BS.pack smallInlineFasta +-- where go (HeaderSize h) (OverlapSize o) (CurrentSize c) = yield (h,o,c) +-- +--smallTest333 = testCase "3/3/3" $ do +-- let res :> r = smallTest 3 3 3 +-- assertEqual "return is null" () r +-- assertEqual "length is 4" 4 (P.length res) +-- assertEqual "!!0" (BioSequenceWindow "Aaa" "" "123" "" (PartialLocation PlusStrand 0 3)) (res!!0) +-- assertEqual "!!1" (BioSequenceWindow "Bbb" "" "456" "" (PartialLocation PlusStrand 0 3)) (res!!1) +-- assertEqual "!!2" (BioSequenceWindow "Bbb" "456" "7" "" (PartialLocation PlusStrand 3 1)) (res!!2) +-- assertEqual "!!3" (BioSequenceWindow "Ccc" "" "890" "" (PartialLocation PlusStrand 0 3)) (res!!3) +-- -- +-- assertEqual "!!0/Fasta" (Fasta "Aaa" "123") (view windowedFasta $ res!!0) +-- assertEqual "!!1/Fasta" (Fasta "Bbb" "456") (view windowedFasta $ res!!1) +-- assertEqual "!!2/Fasta" (Fasta "Bbb" "7" ) (view windowedFasta $ res!!2) +-- assertEqual "!!3/Fasta" (Fasta "Ccc" "890") (view windowedFasta $ res!!3) main :: IO () main = do - gs ← goldenTests - defaultMain $ testGroup "all tests" - [ testGroup "Golden" [gs] - , testGroup "unit tests" [smallTest333] - ] +-- gs ← goldenTests + defaultMain $ testGroup "all tests" +-- [ testGroup "Golden" [gs] +-- [ testGroup "unit tests" [smallTest333] + [ + ]