BiobaseFasta 0.2.0.0 → 0.4.0.1
raw patch · 9 files changed
Files
- Biobase/Fasta/Export.hs +0/−38
- Biobase/Fasta/Streaming.hs +185/−182
- Biobase/Fasta/Strict.hs +107/−0
- Biobase/Fasta/Types.hs +0/−56
- BiobaseFasta.cabal +45/−30
- README.md +2/−1
- changelog.md +6/−0
- src/fastaextract.hs +60/−0
- tests/properties.hs +78/−64
− Biobase/Fasta/Export.hs
@@ -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
Biobase/Fasta/Streaming.hs view
@@ -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+-}
+ Biobase/Fasta/Strict.hs view
@@ -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+
− Biobase/Fasta/Types.hs
@@ -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-
BiobaseFasta.cabal view
@@ -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' <http://hackage.haskell.org/package/fasta> 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
README.md view
@@ -1,4 +1,5 @@-[](https://travis-ci.org/choener/BiobaseFasta)++ # BiobaseFasta
changelog.md view
@@ -1,3 +1,9 @@+0.4.0.1+-------++- streaming-bytestring >= 0.2+- cleanup of older code+ 0.2.0.0 -------
+ src/fastaextract.hs view
@@ -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+
tests/properties.hs view
@@ -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]+ [+ ]