sequence-formats 1.6.3 → 1.6.6.0
raw patch · 10 files changed
+211/−34 lines, 10 files
Files
- Changelog.md +3/−0
- sequence-formats.cabal +5/−1
- src/SequenceFormats/Bed.hs +18/−0
- src/SequenceFormats/Fasta.hs +18/−17
- src/SequenceFormats/Genomic.hs +77/−0
- src/SequenceFormats/Utils.hs +1/−1
- test/SequenceFormats/BedSpec.hs +53/−0
- test/SequenceFormats/FastaSpec.hs +22/−10
- testDat/example.bed +3/−0
- testDat/example.fasta +11/−5
Changelog.md view
@@ -1,3 +1,6 @@+V 1.6.6.0: Added a minimal Bed parser and a bed-filter+V 1.6.5.0: Long chromosome names are now parsed correctly+V 1.6.4: Fixed a bug in Fasta parsing which would skip the first line if the header was minimal. V 1.6.3: Pileup reference-base is read as upper case, even if it's not in the file. V 1.6.2: Made compatible with GHC 9 (thanks to github user sjakobi)
sequence-formats.cabal view
@@ -1,6 +1,6 @@ cabal-version: >=1.10 name: sequence-formats-version: 1.6.3+version: 1.6.6.0 license: GPL-3 license-file: LICENSE maintainer: stephan.schiffels@mac.com@@ -19,6 +19,7 @@ README.md Changelog.md testDat/example.bim+ testDat/example.bed testDat/example.eigenstratgeno testDat/example.fasta testDat/example.freqsum@@ -42,6 +43,8 @@ SequenceFormats.Plink SequenceFormats.Utils SequenceFormats.Pileup+ SequenceFormats.Bed+ SequenceFormats.Genomic hs-source-dirs: src default-language: Haskell2010@@ -67,6 +70,7 @@ hs-source-dirs: test other-modules: SequenceFormats.EigenstratSpec+ SequenceFormats.BedSpec SequenceFormats.FastaSpec SequenceFormats.FreqSumSpec SequenceFormats.RareAlleleHistogramSpec
+ src/SequenceFormats/Bed.hs view
@@ -0,0 +1,18 @@+module SequenceFormats.Bed (BedEntry(..), readBedFile) where++import SequenceFormats.Utils (Chrom(..), readFileProd, consumeProducer)++import Data.Char (isSpace)+import qualified Data.Attoparsec.ByteString.Char8 as A+import Pipes (Producer)+import Pipes.Safe (MonadSafe)++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+ where+ chrom = Chrom <$> A.takeTill isSpace++readBedFile :: (MonadSafe m) => FilePath -> Producer BedEntry m ()+readBedFile bedFile = consumeProducer bedFileParser (readFileProd bedFile)
src/SequenceFormats/Fasta.hs view
@@ -5,24 +5,26 @@ -} module SequenceFormats.Fasta (readNextFastaEntry, loadFastaChrom) where -import SequenceFormats.Utils (Chrom(..))+import SequenceFormats.Utils (Chrom (..)) -import Control.Exception.Base (throwIO, AssertionFailed(..))-import Control.Monad (void)-import Control.Monad.IO.Class (liftIO, MonadIO)-import Control.Monad.Trans.State.Strict (runStateT)+import Control.Exception.Base (AssertionFailed (..),+ throwIO)+import Control.Monad (void)+import Control.Monad.IO.Class (MonadIO, liftIO)+import Control.Monad.Trans.State.Strict (runStateT) import qualified Data.Attoparsec.ByteString.Char8 as A-import qualified Data.ByteString.Char8 as B-import Data.Char (isAlphaNum)-import Lens.Family2 (view)-import Pipes (Producer, next, (>->), runEffect)-import Pipes.Attoparsec (parse)-import qualified Pipes.ByteString as P-import Pipes.Prelude (drain)-import System.IO (Handle, hPutStr, stderr)+import qualified Data.ByteString.Char8 as B+import Data.Char (isSpace)+import Lens.Family2 (view)+import Pipes (Producer, next, runEffect,+ (>->))+import Pipes.Attoparsec (parse)+import qualified Pipes.ByteString as P+import Pipes.Prelude (drain)+import System.IO (Handle, hPutStr, stderr) --- |A function to select out a specific chromosome from a Fasta File. Expects a file handle to the --- file and a chromosome. Note that by Chromosome I simply denote a fasta header line, as is the +-- |A function to select out a specific chromosome from a Fasta File. Expects a file handle to the+-- file and a chromosome. Note that by Chromosome I simply denote a fasta header line, as is the -- case for example for the human reference genome. Returns a Bytestring-Producer over the single sequence followed the specified header (the chromosome). loadFastaChrom :: Handle -> Chrom -> IO (Producer B.ByteString IO ()) loadFastaChrom refFileHandle chrom = do@@ -59,8 +61,7 @@ fastaHeaderLineParser :: A.Parser Chrom fastaHeaderLineParser = do _ <- A.char '>'- chrom <- A.takeWhile isAlphaNum- A.skipSpace+ chrom <- A.takeWhile $ not . isSpace A.skipWhile (\c -> c /= '\n' && c /= '\r') A.endOfLine return . Chrom $ chrom
+ src/SequenceFormats/Genomic.hs view
@@ -0,0 +1,77 @@+module SequenceFormats.Genomic where++import SequenceFormats.Bed (BedEntry(..))+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)++class Genomic a where+ genomicPosition :: a -> (Chrom, Int)++ genomicChrom :: a -> Chrom+ genomicChrom = fst . genomicPosition++ genomicBase :: a -> Int + genomicBase = snd . genomicPosition++instance Genomic EigenstratSnpEntry where+ genomicPosition (EigenstratSnpEntry c p _ _ _ _) = (c, p)++instance Genomic FreqSumEntry where+ genomicPosition (FreqSumEntry c p _ _ _ _ _) = (c, p)++instance Genomic PileupRow where+ genomicPosition (PileupRow c p _ _ _) = (c, p)++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
src/SequenceFormats/Utils.hs view
@@ -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 (\h -> PB.fromHandle h)+readFileProd f = PS.withFile f ReadMode PB.fromHandle word :: A.Parser B.ByteString word = A.takeTill isSpace
+ test/SequenceFormats/BedSpec.hs view
@@ -0,0 +1,53 @@+{-# LANGUAGE OverloadedStrings #-}+module SequenceFormats.BedSpec (spec) where++import Control.Foldl (purely, list)+import Pipes (each)+import qualified Pipes.Prelude as P+import Pipes.Safe (runSafeT)+import SequenceFormats.Bed (BedEntry(..), readBedFile)+import SequenceFormats.Genomic (filterThroughBed)+import SequenceFormats.Eigenstrat (EigenstratSnpEntry(EigenstratSnpEntry))+import SequenceFormats.Utils (Chrom(..))+import Test.Hspec++spec :: Spec+spec = do+ testReadBed+ testFilterThroughBed++mockDatBed :: [BedEntry]+mockDatBed = [+ BedEntry (Chrom "11") 0 100,+ BedEntry (Chrom "11") 200 300,+ BedEntry (Chrom "11") 400 500]+ +testReadBed :: Spec+testReadBed = describe "readBed" $ do+ it "should read the correct eigenstrat file" $ do+ bedEntries <- runSafeT $ purely P.fold list (readBedFile "testDat/example.bed")+ bedEntries `shouldBe` mockDatBed ++mockDatEigenstratSnp :: [EigenstratSnpEntry]+mockDatEigenstratSnp = [+ EigenstratSnpEntry (Chrom "1") 0 0.000000 "rs0000" 'A' 'C',+ EigenstratSnpEntry (Chrom "11") 10 0.001000 "rs1111" 'A' 'G',+ EigenstratSnpEntry (Chrom "11") 100 0.002000 "rs2222" 'A' 'T',+ EigenstratSnpEntry (Chrom "11") 110 0.003000 "rs3333" 'C' 'A',+ EigenstratSnpEntry (Chrom "11") 210 0.004000 "rs4444" 'G' 'A',+ EigenstratSnpEntry (Chrom "11") 310 0.005000 "rs5555" 'T' 'A',+ EigenstratSnpEntry (Chrom "14") 400 0.006000 "rs6666" 'G' 'T']++mockFiltered :: [EigenstratSnpEntry]+mockFiltered = [+ EigenstratSnpEntry (Chrom "11") 10 0.001000 "rs1111" 'A' 'G',+ EigenstratSnpEntry (Chrom "11") 100 0.002000 "rs2222" 'A' 'T',+ EigenstratSnpEntry (Chrom "11") 210 0.004000 "rs4444" 'G' 'A']++testFilterThroughBed :: Spec+testFilterThroughBed = describe "filterThroughBed" $ do+ 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 `shouldBe` mockFiltered
test/SequenceFormats/FastaSpec.hs view
@@ -2,11 +2,11 @@ module SequenceFormats.FastaSpec (spec) where import qualified Data.ByteString.Char8 as BS-import qualified Pipes.Prelude as P-import SequenceFormats.Fasta (loadFastaChrom)-import SequenceFormats.Utils (Chrom(..))-import System.IO (withFile, IOMode(..))-import Test.Hspec+import qualified Pipes.Prelude as P+import SequenceFormats.Fasta (loadFastaChrom)+import SequenceFormats.Utils (Chrom (..))+import System.IO (IOMode (..), withFile)+import Test.Hspec spec :: Spec spec = do@@ -15,9 +15,21 @@ testLoadFastaChrom :: Spec testLoadFastaChrom = around (withFile "testDat/example.fasta" ReadMode) $ describe "loadFastaChrom" $- it "should read the correct fasta string from file" $ \h -> do- fastaProd <- loadFastaChrom h (Chrom "chr3")- fastaString <- BS.unpack . BS.concat <$> P.toListM fastaProd- fastaString `shouldBe` testString+ do+ fastaSpec "short chrom" "chr1"+ "ACGACGACGACGGGGTTTAAAAAGGGTTTCCTCTCTCTCTGGG"++ fastaSpec "short chrom with junk after it" "chr2"+ "ACCAATTTCCCTTTAATATAAGACCCTTTCGGGGAAA"++ fastaSpec "long chrom" "chr4_random"+ "TTGAGAAAGAAATAATAATTTTATATAATTTAAAGGTTATTTAA"++ fastaSpec "long chrom with junk after it" "chr5_random"+ "AAAACTGAAAAAATGTGGTCACCTATGTTAGAACAACAAGGTTT" where- testString = "ACGACGACGACGGGGTTTAAAAAGGGTTTCCTCTCTCTCTGGG"+ fastaSpec desc chr testseq = it (fmtDesc desc) $ \h -> do+ fastaProd <- loadFastaChrom h (Chrom chr)+ fastaString <- BS.unpack . BS.concat <$> P.toListM fastaProd+ fastaString `shouldBe` testseq+ fmtDesc s = "should read the correct fasta string from file (" ++ s ++ ")"
+ testDat/example.bed view
@@ -0,0 +1,3 @@+11 0 100+11 200 300+11 400 500
testDat/example.fasta view
@@ -1,10 +1,16 @@->chr1 bla bla bla+>chr1+ACGACGACGACGGGGTTTAAA+AAGGGTTTCCTCTCTCTCTGGG+>chr2 bla bla bla ACCAATTTCCCTTTAATATAAG ACCCTTTCGGGGAAA->chr2 dsasdfas+>chr3 dsasdfas ACTACTACTACTACTACTACTACGGGGG AAAAAAAAACCCCCCCCCGGGGGGGG->chr3 kflfksja-ACGACGACGACGGGGTTTAAA-AAGGGTTTCCTCTCTCTCTGGG+>chr4_random+TTGAGAAAGAAATAATAATTTT+ATATAATTTAAAGGTTATTTAA+>chr5_random bla bla bla+AAAACTGAAAAAATGTGGTCAC+CTATGTTAGAACAACAAGGTTT