diff --git a/sequence-formats.cabal b/sequence-formats.cabal --- a/sequence-formats.cabal +++ b/sequence-formats.cabal @@ -1,6 +1,6 @@ cabal-version: >=1.10 name: sequence-formats -version: 1.6.6.0 +version: 1.6.6.1 license: GPL-3 license-file: LICENSE maintainer: stephan.schiffels@mac.com @@ -50,19 +50,19 @@ default-language: Haskell2010 build-depends: base >=4.7 && <5, - containers >=0.6.2.1, + containers >=0.6.4.1, errors >=2.3.0, - attoparsec >=0.13.2.4, - pipes >=4.3.14, + attoparsec >=0.14.4, + pipes >=4.3.16, transformers >=0.5.6.2, - bytestring >=0.10.12.0, - lens-family >=2.0.0, - pipes-bytestring >=2.1.6, - foldl >=1.4.10, + bytestring >=0.10.12.1, + lens-family >=2.1.1, + pipes-bytestring >=2.1.7, + foldl >=1.4.12, exceptions >=0.10.4, - pipes-safe >=2.3.2, + pipes-safe >=2.3.3, pipes-attoparsec >=0.5.1.5, - vector >=0.12.1.2 + vector >=0.12.3.1 test-suite sequenceFormatTests type: exitcode-stdio-1.0 @@ -81,15 +81,15 @@ default-language: Haskell2010 build-depends: - base >=4.14.1.0, + base >=4.15.1.0, sequence-formats -any, - foldl >=1.4.10, - pipes >=4.3.14, - pipes-safe >=2.3.2, - tasty >=1.2.3, - vector >=0.12.1.2, + foldl >=1.4.12, + pipes >=4.3.16, + pipes-safe >=2.3.3, + tasty >=1.4.2.1, + vector >=0.12.3.1, transformers >=0.5.6.2, tasty-hunit >=0.10.0.3, - bytestring >=0.10.12.0, - containers >=0.6.2.1, - hspec >=2.7.8 + bytestring >=0.10.12.1, + containers >=0.6.4.1, + hspec >=2.8.5 diff --git a/src/SequenceFormats/Bed.hs b/src/SequenceFormats/Bed.hs --- a/src/SequenceFormats/Bed.hs +++ b/src/SequenceFormats/Bed.hs @@ -1,13 +1,16 @@ -module SequenceFormats.Bed (BedEntry(..), readBedFile) where +module SequenceFormats.Bed where -import SequenceFormats.Utils (Chrom(..), readFileProd, consumeProducer) +import SequenceFormats.Utils (Chrom (..), consumeProducer, + readFileProd) -import Data.Char (isSpace) +import Control.Monad.Trans.Class (lift) import qualified Data.Attoparsec.ByteString.Char8 as A -import Pipes (Producer) -import Pipes.Safe (MonadSafe) +import Data.Char (isSpace) +import Pipes (Producer, next, yield) +import Pipes.Safe (MonadSafe) -data BedEntry = BedEntry Chrom Int Int deriving (Show, Eq) +data BedEntry = BedEntry Chrom Int Int + deriving (Show, Eq) bedFileParser :: A.Parser BedEntry bedFileParser = BedEntry <$> chrom <* A.skipSpace <*> A.decimal <* A.skipSpace <*> A.decimal <* A.endOfLine @@ -16,3 +19,46 @@ readBedFile :: (MonadSafe m) => FilePath -> Producer BedEntry m () readBedFile bedFile = consumeProducer bedFileParser (readFileProd bedFile) + +data IntervalStatus = BedBehind + | BedOn + | BedAhead + +filterThroughBed :: (Monad m) => Producer BedEntry m () -> (a -> (Chrom, Int)) -> Producer a m () -> Producer a m () +filterThroughBed bedProd posFunc gProd = do + b <- lift $ next bedProd + let (bedCurrent, bedRest) = case b of + Left _ -> error "Bed file empty or not readable" + Right r -> r + f' <- lift $ next gProd + let (gCurrent, gRest) = case f' of + Left _ -> error "Genomic stream empty or not readable" + Right r -> r + go bedCurrent gCurrent bedRest gRest + where + go bedCurrent gCurrent bedRest gRest = do + let recurseNextBed = do + b <- lift $ next bedRest + case b of + Left () -> return () + Right (nextBed, bedRest') -> go nextBed gCurrent bedRest' gRest + recurseNextG = do + f' <- lift $ next gRest + case f' of + Left () -> return () + Right (nextG, gRest') -> go bedCurrent nextG bedRest gRest' + case bedCurrent `checkIntervalStatus` gCurrent of + BedBehind -> recurseNextBed + BedAhead -> recurseNextG + BedOn -> do + yield gCurrent + recurseNextG + -- checkIntervalStatus :: BedEntry -> a -> IntervalStatus + checkIntervalStatus (BedEntry bedChrom bedStart bedEnd) g = + case bedChrom `compare` (fst . posFunc) g of + LT -> BedBehind + GT -> BedAhead + EQ -> if bedStart + 1 > (snd . posFunc) g then + BedAhead + else + if bedEnd < (snd . posFunc) g then BedBehind else BedOn diff --git a/src/SequenceFormats/Genomic.hs b/src/SequenceFormats/Genomic.hs --- a/src/SequenceFormats/Genomic.hs +++ b/src/SequenceFormats/Genomic.hs @@ -1,14 +1,13 @@ module SequenceFormats.Genomic where -import SequenceFormats.Bed (BedEntry(..)) +import SequenceFormats.Bed (BedEntry(..), filterThroughBed) import SequenceFormats.Eigenstrat (EigenstratSnpEntry(..)) import SequenceFormats.FreqSum (FreqSumEntry(..)) import SequenceFormats.Pileup (PileupRow(..)) import SequenceFormats.Utils (Chrom) import SequenceFormats.VCF (VCFentry(..)) -import Control.Monad.Trans.Class (lift) -import Pipes (Producer, next, yield) +import Pipes (Producer) class Genomic a where genomicPosition :: a -> (Chrom, Int) @@ -31,47 +30,8 @@ instance Genomic VCFentry where genomicPosition (VCFentry c p _ _ _ _ _ _ _ _) = (c, p) -data IntervalStatus = BedBehind | BedOn | BedAhead - -filterThroughBed :: (Monad m, Genomic e) => Producer BedEntry m () -> Producer e m () -> Producer e m () -filterThroughBed bedProd gProd = do - b <- lift $ next bedProd - let (bedCurrent, bedRest) = case b of - Left _ -> error "Bed file empty or not readable" - Right r -> r - f' <- lift $ next gProd - let (gCurrent, gRest) = case f' of - Left _ -> error "Genomic stream empty or not readable" - Right r -> r - go bedCurrent gCurrent bedRest gRest - where - go bedCurrent gCurrent bedRest gRest = do - let recurseNextBed = do - b <- lift $ next bedRest - case b of - Left () -> return () - Right (nextBed, bedRest') -> go nextBed gCurrent bedRest' gRest - recurseNextG = do - f' <- lift $ next gRest - case f' of - Left () -> return () - Right (nextG, gRest') -> go bedCurrent nextG bedRest gRest' - case bedCurrent `checkIntervalStatus` gCurrent of - BedBehind -> recurseNextBed - BedAhead -> recurseNextG - BedOn -> do - yield gCurrent - recurseNextG - checkIntervalStatus :: (Genomic e) => BedEntry -> e -> IntervalStatus - checkIntervalStatus (BedEntry bedChrom bedStart bedEnd) g = - case bedChrom `compare` genomicChrom g of - LT -> BedBehind - GT -> BedAhead - EQ -> if bedStart + 1 > genomicBase g then - BedAhead - else - if bedEnd < genomicBase g then BedBehind else BedOn - - chromFilter :: (Genomic e) => [Chrom] -> e -> Bool chromFilter exclusionList = (`notElem` exclusionList) . genomicChrom + +genomicFilterThroughBed :: (Monad m, Genomic e) => Producer BedEntry m () -> Producer e m () -> Producer e m () +genomicFilterThroughBed bedProd = filterThroughBed bedProd genomicPosition diff --git a/src/SequenceFormats/Plink.hs b/src/SequenceFormats/Plink.hs --- a/src/SequenceFormats/Plink.hs +++ b/src/SequenceFormats/Plink.hs @@ -13,7 +13,6 @@ import Control.Monad (forM_, void) import Control.Monad.Catch (MonadThrow, throwM) import Control.Monad.IO.Class (MonadIO, liftIO) -import Control.Monad.Trans.Class (lift) import Control.Monad.Trans.State.Strict (runStateT) import qualified Data.Attoparsec.ByteString as AB import qualified Data.Attoparsec.ByteString.Char8 as A @@ -57,8 +56,8 @@ A.skipMany A.space pop <- word ind <- A.skipMany1 A.space >> word - _ <- A.skipMany1 A.space >> A.decimal - _ <- A.skipMany1 A.space >> A.decimal + _ <- A.skipMany1 A.space >> A.decimal :: A.Parser Int + _ <- A.skipMany1 A.space >> A.decimal :: A.Parser Int sex <- A.skipMany1 A.space >> parseSex _ <- A.skipMany1 A.space >> word void A.endOfLine diff --git a/src/SequenceFormats/Utils.hs b/src/SequenceFormats/Utils.hs --- a/src/SequenceFormats/Utils.hs +++ b/src/SequenceFormats/Utils.hs @@ -73,7 +73,7 @@ consumeProducer parser prod = parsed parser prod >>= liftParsingErrors readFileProd :: (PS.MonadSafe m) => FilePath -> Producer B.ByteString m () -readFileProd f = PS.withFile f ReadMode PB.fromHandle +readFileProd f = PS.withFile f ReadMode (\h -> PB.fromHandle h) word :: A.Parser B.ByteString word = A.takeTill isSpace diff --git a/test/SequenceFormats/BedSpec.hs b/test/SequenceFormats/BedSpec.hs --- a/test/SequenceFormats/BedSpec.hs +++ b/test/SequenceFormats/BedSpec.hs @@ -6,7 +6,7 @@ import qualified Pipes.Prelude as P import Pipes.Safe (runSafeT) import SequenceFormats.Bed (BedEntry(..), readBedFile) -import SequenceFormats.Genomic (filterThroughBed) +import SequenceFormats.Genomic (genomicFilterThroughBed) import SequenceFormats.Eigenstrat (EigenstratSnpEntry(EigenstratSnpEntry)) import SequenceFormats.Utils (Chrom(..)) import Test.Hspec @@ -49,5 +49,5 @@ it "should filter correctly" $ do let bedProd = readBedFile "testDat/example.bed" let eigenstratProd = each mockDatEigenstratSnp - filtered <- runSafeT $ purely P.fold list (filterThroughBed bedProd eigenstratProd) + filtered <- runSafeT $ purely P.fold list (genomicFilterThroughBed bedProd eigenstratProd) filtered `shouldBe` mockFiltered