sequence-formats 1.8.0.1 → 1.8.1.0
raw patch · 20 files changed
+337/−300 lines, 20 filesdep +pipes-zlibPVP: major bump suggested
API removals or changes: PVP suggests a major version bump
Dependencies added: pipes-zlib
API changes (from Hackage documentation)
- SequenceFormats.VCF: readVCFfromProd :: MonadThrow m => Producer ByteString m () -> m (VCFheader, Producer VCFentry m ())
+ SequenceFormats.Utils: readFileProdCheckCompress :: MonadSafe m => FilePath -> Producer ByteString m ()
Files
- Changelog.md +1/−0
- sequence-formats.cabal +60/−91
- src/SequenceFormats/Bed.hs +1/−1
- src/SequenceFormats/Eigenstrat.hs +8/−4
- src/SequenceFormats/FreqSum.hs +28/−27
- src/SequenceFormats/Genomic.hs +8/−8
- src/SequenceFormats/Pileup.hs +15/−14
- src/SequenceFormats/Plink.hs +4/−3
- src/SequenceFormats/RareAlleleHistogram.hs +29/−28
- src/SequenceFormats/Utils.hs +23/−16
- src/SequenceFormats/VCF.hs +36/−34
- test/SequenceFormats/BedSpec.hs +11/−11
- test/SequenceFormats/EigenstratSpec.hs +30/−14
- test/SequenceFormats/FreqSumSpec.hs +13/−13
- test/SequenceFormats/PileupSpec.hs +8/−7
- test/SequenceFormats/PlinkSpec.hs +25/−6
- test/SequenceFormats/RareAlleleHistogramSpec.hs +5/−3
- test/SequenceFormats/UtilsSpec.hs +4/−4
- test/SequenceFormats/VCFSpec.hs +28/−9
- testDat/example.bim +0/−7
Changelog.md view
@@ -1,5 +1,6 @@ # Changelog +- V 1.8.1.0: Added gzip-support (read-only for now) for Plink (bed and bim files) and VCF. - V 1.8.0.1: Allow reading arbitrary letters as reference base in Pileup format. Before only ACTG, N and M were allowed. Now all letters are allowed, as we found out that reference fasta files can contain more letters than just those and we don't want to immediately stop parsing in that case.
sequence-formats.cabal view
@@ -1,95 +1,64 @@-cabal-version: >=1.10-name: sequence-formats-version: 1.8.0.1-license: GPL-3-license-file: LICENSE-maintainer: stephan.schiffels@mac.com-author: Stephan Schiffels-homepage: https://github.com/stschiff/sequence-formats-bug-reports: https://github.com/stschiff/sequence-formats/issues-synopsis:- A package with basic parsing utilities for several Bioinformatic data formats.--description:- Contains utilities to parse and write Eigenstrat, Fasta, FreqSum, VCF, Plink and other file formats used in population genetics analyses.+name: sequence-formats+version: 1.8.1.0+synopsis: A package with basic parsing utilities for several Bioinformatic data formats.+description: Contains utilities to parse and write Eigenstrat, Fasta, FreqSum, VCF, Plink and other file formats used in population genetics analyses.+license: GPL-3+license-file: LICENSE+author: Stephan Schiffels+maintainer: stephan.schiffels@mac.com+category: Bioinformatics+build-type: Simple+cabal-version: >=1.10+Homepage: https://github.com/stschiff/sequence-formats+Bug-Reports: https://github.com/stschiff/sequence-formats/issues -category: Bioinformatics-build-type: Simple-extra-source-files:- README.md- Changelog.md- testDat/example.bim- testDat/example.bed- testDat/example.eigenstratgeno- testDat/example.fasta- testDat/example.freqsum- testDat/example.histogram.txt- testDat/example.ind- testDat/example.snp- testDat/example.vcf- testDat/example.pileup- testDat/example.fam- testDat/example.plink.bed- testDat/example.plink.fam- testDat/example.plink.bim+extra-source-files: README.md,+ Changelog.md,+ testDat/example.bed+ testDat/example.eigenstratgeno,+ testDat/example.fasta,+ testDat/example.freqsum,+ testDat/example.histogram.txt,+ testDat/example.ind,+ testDat/example.snp,+ testDat/example.vcf,+ testDat/example.pileup+ testDat/example.fam+ testDat/example.plink.bed+ testDat/example.plink.fam+ testDat/example.plink.bim library- exposed-modules:- SequenceFormats.RareAlleleHistogram- SequenceFormats.FreqSum- SequenceFormats.Fasta- SequenceFormats.VCF- SequenceFormats.Eigenstrat- SequenceFormats.Plink- SequenceFormats.Utils- SequenceFormats.Pileup- SequenceFormats.Bed- SequenceFormats.Genomic-- hs-source-dirs: src- default-language: Haskell2010- build-depends:- base >=4.7 && <5,- containers,- errors,- attoparsec,- pipes,- transformers,- bytestring,- lens-family,- pipes-bytestring,- foldl,- exceptions,- pipes-safe,- pipes-attoparsec,- vector--test-suite sequenceFormatTests- type: exitcode-stdio-1.0- main-is: Spec.hs- hs-source-dirs: test- other-modules:- SequenceFormats.EigenstratSpec- SequenceFormats.BedSpec- SequenceFormats.FastaSpec- SequenceFormats.FreqSumSpec- SequenceFormats.RareAlleleHistogramSpec- SequenceFormats.UtilsSpec- SequenceFormats.PlinkSpec- SequenceFormats.VCFSpec- SequenceFormats.PileupSpec+ exposed-modules: SequenceFormats.RareAlleleHistogram,+ SequenceFormats.FreqSum,+ SequenceFormats.Fasta,+ SequenceFormats.VCF,+ SequenceFormats.Eigenstrat,+ SequenceFormats.Plink,+ SequenceFormats.Utils,+ SequenceFormats.Pileup,+ SequenceFormats.Bed,+ SequenceFormats.Genomic+ hs-source-dirs: src+ build-depends: base >= 4.7 && < 5, containers, errors, attoparsec, pipes,+ transformers, bytestring, lens-family,+ pipes-bytestring, foldl, exceptions, pipes-safe,+ pipes-attoparsec, vector, pipes-zlib+ default-language: Haskell2010 - default-language: Haskell2010- build-depends:- base,- sequence-formats,- foldl,- pipes,- pipes-safe,- tasty,- vector,- transformers,- tasty-hunit,- bytestring,- containers,- hspec+Test-Suite sequenceFormatTests+ type: exitcode-stdio-1.0+ main-is: Spec.hs+ hs-source-dirs: test+ build-depends: base, sequence-formats, foldl, pipes, pipes-safe, tasty, vector,+ transformers, tasty-hunit, bytestring, containers, hspec+ other-modules: SequenceFormats.EigenstratSpec,+ SequenceFormats.BedSpec,+ SequenceFormats.FastaSpec,+ SequenceFormats.FreqSumSpec,+ SequenceFormats.RareAlleleHistogramSpec,+ SequenceFormats.UtilsSpec,+ SequenceFormats.PlinkSpec,+ SequenceFormats.VCFSpec,+ SequenceFormats.PileupSpec+ default-language: Haskell2010
src/SequenceFormats/Bed.hs view
@@ -45,7 +45,7 @@ recurseNextG = do f' <- lift $ next gRest case f' of- Left () -> return ()+ Left () -> return () Right (nextG, gRest') -> go bedCurrent nextG bedRest gRest' case bedCurrent `checkIntervalStatus` gCurrent of BedBehind -> recurseNextBed
src/SequenceFormats/Eigenstrat.hs view
@@ -13,7 +13,8 @@ import SequenceFormats.Utils (Chrom (..), SeqFormatException (..), consumeProducer,- readFileProd, word)+ readFileProdCheckCompress,+ word) import Control.Applicative ((<|>)) import Control.Exception (throw)@@ -121,8 +122,12 @@ -- |Function to read a Snp File from a file. Returns a Pipes-Producer over the EigenstratSnpEntries. readEigenstratSnpFile :: (MonadSafe m) => FilePath -> Producer EigenstratSnpEntry m ()-readEigenstratSnpFile = consumeProducer eigenstratSnpParser . readFileProd+readEigenstratSnpFile = consumeProducer eigenstratSnpParser . readFileProdCheckCompress +-- |Function to read a Geno File from a file. Returns a Pipes-Producer over the GenoLines.+readEigenstratGenoFile :: (MonadSafe m) => FilePath -> Producer GenoLine m ()+readEigenstratGenoFile = consumeProducer eigenstratGenoParser . readFileProdCheckCompress+ -- |Function to read a full Eigenstrat database from files. Returns a pair of the Eigenstrat Individual Entries, and a joint Producer over the snp entries and the genotypes. readEigenstrat :: (MonadSafe m) => FilePath -- ^The Genotype file -> FilePath -- ^The Snp File@@ -131,8 +136,7 @@ readEigenstrat genoFile snpFile indFile = do indEntries <- readEigenstratInd indFile let snpProd = readEigenstratSnpFile snpFile- genoProd = consumeProducer eigenstratGenoParser (readFileProd genoFile) >->- validateEigenstratEntries (length indEntries)+ genoProd = readEigenstratGenoFile genoFile >-> validateEigenstratEntries (length indEntries) return (indEntries, P.zip snpProd genoProd) validateEigenstratEntries :: (MonadThrow m) => Int -> Pipe GenoLine GenoLine m ()
src/SequenceFormats/FreqSum.hs view
@@ -4,30 +4,31 @@ <https://rarecoal-docs.readthedocs.io/en/latest/rarecoal-tools.html#vcf2freqsum> -} -module SequenceFormats.FreqSum (readFreqSumStdIn, readFreqSumFile, FreqSumEntry(..), +module SequenceFormats.FreqSum (readFreqSumStdIn, readFreqSumFile, FreqSumEntry(..), FreqSumHeader(..), printFreqSumStdOut, printFreqSumFile, freqSumEntryToText) where -import SequenceFormats.Utils (consumeProducer, Chrom(..), readFileProd)+import SequenceFormats.Utils (Chrom (..), consumeProducer,+ readFileProd) -import Control.Applicative ((<|>))-import Control.Monad.Catch (MonadThrow, throwM)-import Control.Monad.IO.Class (MonadIO, liftIO)-import Control.Monad.Trans.State.Strict (runStateT)+import Control.Applicative ((<|>))+import Control.Monad.Catch (MonadThrow, throwM)+import Control.Monad.IO.Class (MonadIO, liftIO)+import Control.Monad.Trans.State.Strict (runStateT) import qualified Data.Attoparsec.ByteString.Char8 as A-import Data.Char (isAlphaNum, isSpace)-import qualified Data.ByteString.Char8 as B-import Pipes (Producer, (>->), Consumer)-import Pipes.Attoparsec (parse, ParsingError(..))-import qualified Pipes.Prelude as P-import Pipes.Safe (MonadSafe)-import Pipes.Safe.Prelude (withFile)-import qualified Pipes.ByteString as PB-import Prelude hiding (putStr)-import System.IO (IOMode(..))+import qualified Data.ByteString.Char8 as B+import Data.Char (isAlphaNum, isSpace)+import Pipes (Consumer, Producer, (>->))+import Pipes.Attoparsec (ParsingError (..), parse)+import qualified Pipes.ByteString as PB+import qualified Pipes.Prelude as P+import Pipes.Safe (MonadSafe)+import Pipes.Safe.Prelude (withFile)+import Prelude hiding (putStr)+import System.IO (IOMode (..)) -- |A Datatype representing the Header data FreqSumHeader = FreqSumHeader {- fshNames :: [String], -- ^A list of individual or group names+ fshNames :: [String], -- ^A list of individual or group names fshCounts :: [Int] -- ^A list of haplotype counts per individual/group. } deriving (Eq, Show) @@ -39,13 +40,13 @@ -- |A Datatype to denote a single freqSum line data FreqSumEntry = FreqSumEntry {- fsChrom :: Chrom, -- ^The chromosome of the site- fsPos :: Int, -- ^The position of the site- fsSnpId :: Maybe B.ByteString, -- ^An optional parameter to take the snpId. This is not parsed from or printed to freqSum format but is used in internal conversions from Eigenstrat.+ fsChrom :: Chrom, -- ^The chromosome of the site+ fsPos :: Int, -- ^The position of the site+ fsSnpId :: Maybe B.ByteString, -- ^An optional parameter to take the snpId. This is not parsed from or printed to freqSum format but is used in internal conversions from Eigenstrat. fsGeneticPos :: Maybe Double, -- ^An optional parameter to take the genetic pos. This is not parsed from or printed to freqSum format but is used in internal conversions from Eigenstrat.- fsRef :: Char, -- ^The reference allele- fsAlt :: Char, -- ^The alternative allele- fsCounts :: [Maybe Int] -- ^A list of allele counts in each group. Nothing denotes missing data.+ fsRef :: Char, -- ^The reference allele+ fsAlt :: Char, -- ^The alternative allele+ fsCounts :: [Maybe Int] -- ^A list of allele counts in each group. Nothing denotes missing data. } deriving (Eq, Show) -- |This function converts a single freqSum entry to a printable freqSum line.@@ -53,8 +54,8 @@ freqSumEntryToText (FreqSumEntry chrom pos _ _ ref alt maybeCounts) = B.intercalate "\t" [unChrom chrom, B.pack (show pos), B.singleton ref, B.singleton alt, countStr] <> "\n" where- countStr = B.intercalate "\t" . map (B.pack . show . convertToNum) $ maybeCounts - convertToNum Nothing = -1+ countStr = B.intercalate "\t" . map (B.pack . show . convertToNum) $ maybeCounts+ convertToNum Nothing = -1 convertToNum (Just a) = a readFreqSumProd :: (MonadThrow m) =>@@ -62,8 +63,8 @@ readFreqSumProd prod = do (res, rest) <- runStateT (parse parseFreqSumHeader) prod header <- case res of- Nothing -> throwM $ ParsingError [] "freqSum file exhausted"- Just (Left e) -> throwM e+ Nothing -> throwM $ ParsingError [] "freqSum file exhausted"+ Just (Left e) -> throwM e Just (Right h) -> return h return (header, consumeProducer parseFreqSumEntry rest)
src/SequenceFormats/Genomic.hs view
@@ -1,13 +1,13 @@ module SequenceFormats.Genomic where -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 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 Pipes (Producer)+import Pipes (Producer) class Genomic a where genomicPosition :: a -> (Chrom, Int)@@ -15,7 +15,7 @@ genomicChrom :: a -> Chrom genomicChrom = fst . genomicPosition - genomicBase :: a -> Int + genomicBase :: a -> Int genomicBase = snd . genomicPosition instance Genomic EigenstratSnpEntry where
src/SequenceFormats/Pileup.hs view
@@ -1,16 +1,17 @@ {-# LANGUAGE OverloadedStrings #-} module SequenceFormats.Pileup (readPileupFromStdIn, readPileupFromFile, PileupRow(..), Strand(..)) where -import Control.Monad.Catch (MonadThrow)-import Control.Monad.IO.Class (MonadIO)+import Control.Monad.Catch (MonadThrow)+import Control.Monad.IO.Class (MonadIO) import qualified Data.Attoparsec.ByteString.Char8 as A-import qualified Data.ByteString.Char8 as B-import Data.Char (toUpper)-import Pipes (Producer)-import qualified Pipes.ByteString as PB-import Pipes.Safe (MonadSafe)+import qualified Data.ByteString.Char8 as B+import Data.Char (toUpper)+import Pipes (Producer)+import qualified Pipes.ByteString as PB+import Pipes.Safe (MonadSafe) -import SequenceFormats.Utils (Chrom(..), word, readFileProd, consumeProducer)+import SequenceFormats.Utils (Chrom (..), consumeProducer,+ readFileProd, word) -- |A datatype to represent the strand orientation of a single base. data Strand = ForwardStrand | ReverseStrand deriving (Eq, Show)@@ -19,10 +20,10 @@ -- The constructor arguments are: Chromosome, Position, Refererence Allelele, -- Pileup String per individual data PileupRow = PileupRow {- pileupChrom :: Chrom, -- ^The chromosome- pileupPos :: Int, -- ^The position- pileupRef :: Char, -- ^The reference base- pileupBases :: [String], -- ^The base string+ pileupChrom :: Chrom, -- ^The chromosome+ pileupPos :: Int, -- ^The position+ pileupRef :: Char, -- ^The reference base+ pileupBases :: [String], -- ^The base string pileupStrandInfo :: [[Strand]] } deriving (Eq, Show) @@ -42,7 +43,7 @@ pos <- A.decimal _ <- A.space refA <- toUpper <$> A.letter_ascii- -- we used to parse only ACTGN and M and the small-letter versions of those, + -- we used to parse only ACTGN and M and the small-letter versions of those, -- but it turns out that there can be more letter from the IUPAC alphabet that can occur in -- fasta reference files, so here we're parsing them all in principle. -- Also, for some reason, there is an M in the human reference at@@ -72,7 +73,7 @@ | x `elem` ("ACTGN" :: String) = (x, ForwardStrand) : go xs | x `elem` ("actgn" :: String) = (toUpper x, ReverseStrand) : go xs | x `elem` ("$*#<>" :: String) = go xs- | x == '^' = go (drop 1 xs) -- skip the next character, which is the mapping quality + | x == '^' = go (drop 1 xs) -- skip the next character, which is the mapping quality | x == '+' || x == '-' = -- insertions or deletions, followed by a decimal number case reads xs of [(num, rest)] -> go (drop num rest)
src/SequenceFormats/Plink.hs view
@@ -17,7 +17,8 @@ GenoEntry (..), GenoLine, Sex (..)) import SequenceFormats.Utils (Chrom (..), consumeProducer,- readFileProd, word)+ readFileProdCheckCompress,+ word) import Control.Applicative ((<|>)) import Control.Monad (forM_, void)@@ -137,7 +138,7 @@ -- |A function to read a bed file from a file. Returns a Producer over all lines. readPlinkBedFile :: (MonadSafe m) => FilePath -> Int -> m (Producer GenoLine m ())-readPlinkBedFile file nrInds = readPlinkBedProd nrInds (readFileProd file)+readPlinkBedFile file nrInds = readPlinkBedProd nrInds . readFileProdCheckCompress $ file -- |Function to read a Bim File from StdIn. Returns a Pipes-Producer over the EigenstratSnpEntries. readBimStdIn :: (MonadThrow m, MonadIO m) => Producer EigenstratSnpEntry m ()@@ -145,7 +146,7 @@ -- |Function to read a Bim File from a file. Returns a Pipes-Producer over the EigenstratSnpEntries. readBimFile :: (MonadSafe m) => FilePath -> Producer EigenstratSnpEntry m ()-readBimFile = consumeProducer bimParser . readFileProd+readBimFile = consumeProducer bimParser . readFileProdCheckCompress -- |Function to read a Plink fam file. Returns the Eigenstrat Individual Entries as list. readFamFile :: (MonadIO m) => FilePath -> m [PlinkFamEntry]
src/SequenceFormats/RareAlleleHistogram.hs view
@@ -7,38 +7,39 @@ module SequenceFormats.RareAlleleHistogram (RareAlleleHistogram(..), readHistogramFromHandle, SitePattern, readHistogram, writeHistogramStdOut, writeHistogramFile, showSitePattern) where -import SequenceFormats.Utils (SeqFormatException(..))+import SequenceFormats.Utils (SeqFormatException (..)) -import Control.Applicative (optional)-import Control.Error (assertErr)-import Control.Exception (throw)-import Control.Monad.IO.Class (MonadIO, liftIO)-import Control.Monad.Trans.State.Strict (evalStateT)+import Control.Applicative (optional)+import Control.Error (assertErr)+import Control.Exception (throw)+import Control.Monad.IO.Class (MonadIO, liftIO)+import Control.Monad.Trans.State.Strict (evalStateT) import qualified Data.Attoparsec.ByteString.Char8 as A-import Data.Char (isAlphaNum)-import Data.Int (Int64)-import Data.List (intercalate, sortBy)-import qualified Data.Map.Strict as Map-import qualified Data.ByteString.Char8 as B-import Pipes.Attoparsec (parse)-import qualified Pipes.ByteString as PB-import System.IO (Handle, IOMode(..), withFile)+import qualified Data.ByteString.Char8 as B+import Data.Char (isAlphaNum)+import Data.Int (Int64)+import Data.List (intercalate, sortBy)+import qualified Data.Map.Strict as Map+import Pipes.Attoparsec (parse)+import qualified Pipes.ByteString as PB+import System.IO (Handle, IOMode (..),+ withFile) -- |A datatype to represent an Allele Sharing Histogram: data RareAlleleHistogram = RareAlleleHistogram {- raNames :: [String], -- ^A list of branch names- raNVec :: [Int], -- ^A list of haploid sample sizes.- raMinAf :: Int, -- ^The minimum allele count- raMaxAf :: Int, -- ^The maximum allele count- raConditionOn :: [Int], -- ^A list of branch indices that were used to condition the allele + raNames :: [String], -- ^A list of branch names+ raNVec :: [Int], -- ^A list of haploid sample sizes.+ raMinAf :: Int, -- ^The minimum allele count+ raMaxAf :: Int, -- ^The maximum allele count+ raConditionOn :: [Int], -- ^A list of branch indices that were used to condition the allele --sharing pattern- raExcludePatterns :: [SitePattern], -- ^A list of patterns that are excluded.- raTotalNrSites :: Int64, -- ^The total number of non-missing sites in the genome.- raCounts :: Map.Map SitePattern Int64, -- ^The actual data, a dictionary from allele sharing patterns to observed numbers.+ raExcludePatterns :: [SitePattern], -- ^A list of patterns that are excluded.+ raTotalNrSites :: Int64, -- ^The total number of non-missing sites in the genome.+ raCounts :: Map.Map SitePattern Int64, -- ^The actual data, a dictionary from allele sharing patterns to observed numbers. raJackknifeEstimates :: Maybe (Map.Map SitePattern (Double, Double)) -- ^An optional dictionary that contains Jackknife estimates and standard deviations for each pattern frequency. } deriving (Eq, Show) --- |A simple type synonym for the SitePattern, represented as a list of Integers that represents +-- |A simple type synonym for the SitePattern, represented as a list of Integers that represents -- each pattern across the branches. type SitePattern = [Int] @@ -46,8 +47,8 @@ showSitePattern :: SitePattern -> String showSitePattern = intercalate "," . map show --- |Function to convert a Rare Allele Histogram to text. Returns an error if attempting to print a --- histogram with non-standard settings. Many settings, such as minAf>1, are only meant for +-- |Function to convert a Rare Allele Histogram to text. Returns an error if attempting to print a+-- histogram with non-standard settings. Many settings, such as minAf>1, are only meant for -- in-memory representations, but are not compatible with the file format itself. showHistogram :: RareAlleleHistogram -> Either String B.ByteString showHistogram hist = do@@ -76,16 +77,16 @@ writeHistogramStdOut :: (MonadIO m) => RareAlleleHistogram -> m () writeHistogramStdOut hist = case showHistogram hist of- Left err -> throw (SeqFormatException err)+ Left err -> throw (SeqFormatException err) Right outStr -> liftIO $ B.putStrLn outStr -- |Write a histogram to a file writeHistogramFile :: (MonadIO m) => FilePath -> RareAlleleHistogram -> m () writeHistogramFile outF hist = case showHistogram hist of- Left err -> throw (SeqFormatException err)+ Left err -> throw (SeqFormatException err) Right outStr -> liftIO $ B.writeFile outF outStr- + -- |Read a histogram from a FilePath readHistogram :: (MonadIO m) => FilePath -> m RareAlleleHistogram readHistogram path = liftIO $ withFile path ReadMode readHistogramFromHandle
src/SequenceFormats/Utils.hs view
@@ -3,23 +3,25 @@ -- |This module contains helper functions for file parsing. module SequenceFormats.Utils (liftParsingErrors,- consumeProducer, readFileProd,+ consumeProducer, readFileProd, readFileProdCheckCompress, SeqFormatException(..), Chrom(..), word) where -import Control.Error (readErr)-import Control.Exception (Exception, throw)-import Control.Monad.Catch (MonadThrow, throwM)-import Control.Monad.Trans.Class (lift)+import Control.Error (readErr)+import Control.Exception (Exception, throw)+import Control.Monad.Catch (MonadThrow, throwM)+import Control.Monad.Trans.Class (lift) import qualified Data.Attoparsec.ByteString.Char8 as A-import qualified Data.ByteString.Char8 as B-import Data.Char (isSpace)-import Pipes (Producer, next)-import Pipes.Attoparsec (ParsingError(..), parsed)-import qualified Pipes.ByteString as PB-import qualified Pipes.Safe as PS-import qualified Pipes.Safe.Prelude as PS-import System.IO (IOMode(..))+import qualified Data.ByteString.Char8 as B+import Data.Char (isSpace)+import Data.List (isSuffixOf)+import Pipes (Producer, next)+import Pipes.Attoparsec (ParsingError (..), parsed)+import qualified Pipes.ByteString as PB+import Pipes.GZip (decompress)+import qualified Pipes.Safe as PS+import qualified Pipes.Safe.Prelude as PS+import System.IO (IOMode (..)) -- |An exception type for parsing BioInformatic file formats. data SeqFormatException = SeqFormatException String@@ -36,11 +38,11 @@ -- |Ord instance for Chrom instance Ord Chrom where- compare (Chrom c1) (Chrom c2) = + compare (Chrom c1) (Chrom c2) = let (c1NoChr, c2NoChr) = (removeChr c1, removeChr c2) (c1XYMTconvert, c2XYMTconvert) = (convertXYMT c1NoChr, convertXYMT c2NoChr) in case (,) <$> readChrom c1XYMTconvert <*> readChrom c2XYMTconvert of- Left e -> throw e+ Left e -> throw e Right (cn1, cn2) -> cn1 `compare` cn2 where removeChr :: B.ByteString -> B.ByteString@@ -54,7 +56,7 @@ readChrom :: B.ByteString -> Either SeqFormatException Int readChrom c = readErr (SeqFormatException $ "cannot parse chromosome " ++ B.unpack c) . B.unpack $ c --- |A function to help with reporting parsing errors to stderr. Returns a clean Producer over the +-- |A function to help with reporting parsing errors to stderr. Returns a clean Producer over the -- parsed datatype. liftParsingErrors :: (MonadThrow m) => Either (ParsingError, Producer B.ByteString m r) () -> Producer a m ()@@ -74,6 +76,11 @@ readFileProd :: (PS.MonadSafe m) => FilePath -> Producer B.ByteString m () readFileProd f = PS.withFile f ReadMode (\h -> PB.fromHandle h)++readFileProdCheckCompress :: (PS.MonadSafe m) => FilePath -> Producer B.ByteString m ()+readFileProdCheckCompress f =+ let decompressFunc = if ".gz" `isSuffixOf` f then decompress else id+ in decompressFunc $ PS.withFile f ReadMode (\h -> PB.fromHandle h) word :: A.Parser B.ByteString word = A.takeTill isSpace
src/SequenceFormats/VCF.hs view
@@ -8,48 +8,50 @@ VCFentry(..), readVCFfromStdIn, readVCFfromFile,- readVCFfromProd, getGenotypes, getDosages, isTransversionSnp, vcfToFreqSumEntry, isBiallelicSnp) where -import SequenceFormats.Utils (consumeProducer, Chrom(..),- readFileProd, SeqFormatException(..), word)-import SequenceFormats.FreqSum (FreqSumEntry(..))+import SequenceFormats.FreqSum (FreqSumEntry (..))+import SequenceFormats.Utils (Chrom (..),+ SeqFormatException (..),+ consumeProducer,+ readFileProdCheckCompress,+ word) -import Control.Applicative ((<|>))-import Control.Error (headErr, assertErr)-import Control.Monad (void)-import Control.Monad.Catch (MonadThrow, throwM)-import Control.Monad.Trans.State.Strict (runStateT)-import Control.Monad.IO.Class (MonadIO)+import Control.Applicative ((<|>))+import Control.Error (assertErr, headErr)+import Control.Monad (void)+import Control.Monad.Catch (MonadThrow, throwM)+import Control.Monad.IO.Class (MonadIO)+import Control.Monad.Trans.State.Strict (runStateT) import qualified Data.Attoparsec.ByteString.Char8 as A-import Data.Char (isSpace)-import qualified Data.ByteString.Char8 as B-import Pipes (Producer)-import Pipes.Attoparsec (parse)-import Pipes.Safe (MonadSafe)-import qualified Pipes.ByteString as PB+import qualified Data.ByteString.Char8 as B+import Data.Char (isSpace)+import Pipes (Producer)+import Pipes.Attoparsec (parse)+import qualified Pipes.ByteString as PB+import Pipes.Safe (MonadSafe) -- |A datatype to represent the VCF Header. Most comments are simply parsed as entire lines, but the very last comment line, containing the sample names, is separated out data VCFheader = VCFheader { vcfHeaderComments :: [String], -- ^A list of containing all comments starting with a single '#'- vcfSampleNames :: [String] -- ^The list of sample names parsed from the last comment line + vcfSampleNames :: [String] -- ^The list of sample names parsed from the last comment line -- starting with '##' } deriving (Show) -- |A Datatype representing a single VCF entry. data VCFentry = VCFentry {- vcfChrom :: Chrom, -- ^The chromosome- vcfPos :: Int, -- ^The position- vcfId :: Maybe B.ByteString, -- ^The SNP ID if non-missing- vcfRef :: B.ByteString, -- ^ The reference allele (supports also multi-character alleles for Indels)- vcfAlt :: [B.ByteString], -- ^The alternative alleles, each one possible of multiple characters - vcfQual :: Maybe Double, -- ^The quality value- vcfFilter :: Maybe B.ByteString, -- ^The Filter value, if non-missing.- vcfInfo :: [B.ByteString], -- ^A list of Info fields+ vcfChrom :: Chrom, -- ^The chromosome+ vcfPos :: Int, -- ^The position+ vcfId :: Maybe B.ByteString, -- ^The SNP ID if non-missing+ vcfRef :: B.ByteString, -- ^ The reference allele (supports also multi-character alleles for Indels)+ vcfAlt :: [B.ByteString], -- ^The alternative alleles, each one possible of multiple characters+ vcfQual :: Maybe Double, -- ^The quality value+ vcfFilter :: Maybe B.ByteString, -- ^The Filter value, if non-missing.+ vcfInfo :: [B.ByteString], -- ^A list of Info fields vcfFormatString :: [B.ByteString], -- ^A list of format tags vcfGenotypeInfo :: [[B.ByteString]] -- ^A list of format fields for each sample. } deriving (Show, Eq)@@ -60,8 +62,8 @@ readVCFfromProd prod = do (res, rest) <- runStateT (parse vcfHeaderParser) prod header <- case res of- Nothing -> throwM $ SeqFormatException "freqSum file exhausted"- Just (Left e) -> throwM (SeqFormatException (show e))+ Nothing -> throwM $ SeqFormatException "freqSum file exhausted"+ Just (Left e) -> throwM (SeqFormatException (show e)) Just (Right h) -> return h return (header, consumeProducer vcfEntryParser rest) @@ -71,7 +73,7 @@ -- |Reading a VCF from a file. Returns a VCFHeader and a Producer over VCFentries. readVCFfromFile :: (MonadSafe m) => FilePath -> m (VCFheader, Producer VCFentry m ())-readVCFfromFile = readVCFfromProd . readFileProd+readVCFfromFile = readVCFfromProd . readFileProdCheckCompress vcfHeaderParser :: A.Parser VCFheader vcfHeaderParser = VCFheader <$> A.many1' doubleCommentLine <*> singleCommentLine@@ -92,8 +94,8 @@ vcfEntryParser = vcfEntryParserFull <|> vcfEntryParserTruncated where vcfEntryParserFull = VCFentry <$> (Chrom <$> word) <* sp <*> A.decimal <* sp <*> parseId <*- sp <*> word <* sp <*> parseAlternativeAlleles <* sp <*> parseQual <* sp <*> parseFilter <* - sp <*> parseInfoFields <* sp <*> parseFormatStrings <* sp <*> parseGenotypeInfos <* + sp <*> word <* sp <*> parseAlternativeAlleles <* sp <*> parseQual <* sp <*> parseFilter <*+ sp <*> parseInfoFields <* sp <*> parseFormatStrings <* sp <*> parseGenotypeInfos <* A.endOfLine vcfEntryParserTruncated = VCFentry <$> (Chrom <$> word) <* sp <*> A.decimal <* sp <*> parseId <* sp <*> word <* sp <*> parseAlternativeAlleles <* sp <*> parseQual <* sp <*> parseFilter <*@@ -111,7 +113,7 @@ parseFormatString = A.takeTill (\c -> c == ':' || isSpace c) parseGenotypeInfos = parseGenotype `A.sepBy1` sp parseGenotype = parseGenoField `A.sepBy1` A.char ':'- parseGenoField = A.takeTill (\c -> c == ':' || isSpace c) + parseGenoField = A.takeTill (\c -> c == ':' || isSpace c) -- |returns True if the SNP is biallelic. isBiallelicSnp :: B.ByteString -> [B.ByteString] -> Bool@@ -120,14 +122,14 @@ validRef = (ref `elem` ["A", "C", "G", "T"]) validAlt = case alt of [alt'] -> alt' `elem` ["A", "C", "G", "T"]- _ -> False+ _ -> False -- |returns True if the SNp is a biallelic Transversion SNP (i.e. one of G/T, G/C, A/T, A/C) isTransversionSnp :: B.ByteString -> [B.ByteString] -> Bool isTransversionSnp ref alt = case alt of [alt'] -> isBiallelicSnp ref alt && (not $ isTransition ref alt')- _ -> False+ _ -> False where isTransition r a = ((r == "A") && (a == "G")) || ((r == "G") && (a == "A")) || ((r == "C") && (a == "T")) || ((r == "T") && (a == "C"))@@ -163,4 +165,4 @@ assertErr "Invalid Reference Allele" $ ref `elem` ['A', 'C', 'T', 'G', 'N'] assertErr "Invalid Alternative Allele" $ alt `elem` ['A', 'C', 'T', 'G', '.'] return $ FreqSumEntry (vcfChrom vcfEntry) (vcfPos vcfEntry) (vcfId vcfEntry) Nothing ref alt dosages- +
test/SequenceFormats/BedSpec.hs view
@@ -1,15 +1,15 @@ {-# 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 (genomicFilterThroughBed)-import SequenceFormats.Eigenstrat (EigenstratSnpEntry(EigenstratSnpEntry))-import SequenceFormats.Utils (Chrom(..))-import Test.Hspec+import Control.Foldl (list, purely)+import Pipes (each)+import qualified Pipes.Prelude as P+import Pipes.Safe (runSafeT)+import SequenceFormats.Bed (BedEntry (..), readBedFile)+import SequenceFormats.Eigenstrat (EigenstratSnpEntry (EigenstratSnpEntry))+import SequenceFormats.Genomic (genomicFilterThroughBed)+import SequenceFormats.Utils (Chrom (..))+import Test.Hspec spec :: Spec spec = do@@ -21,12 +21,12 @@ 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 + bedEntries `shouldBe` mockDatBed mockDatEigenstratSnp :: [EigenstratSnpEntry] mockDatEigenstratSnp = [
test/SequenceFormats/EigenstratSpec.hs view
@@ -1,20 +1,23 @@ {-# LANGUAGE OverloadedStrings #-} module SequenceFormats.EigenstratSpec (spec) where -import Control.Foldl (purely, list)-import Control.Monad.IO.Class (liftIO)-import Data.Vector (fromList)-import Pipes (each, runEffect, (>->))-import qualified Pipes.Prelude as P-import Pipes.Safe (runSafeT)-import SequenceFormats.Eigenstrat (readEigenstrat, writeEigenstrat,- EigenstratSnpEntry(..), EigenstratIndEntry(..), GenoLine, Sex(..), GenoEntry(..))-import SequenceFormats.Utils (Chrom(..))-import Test.Hspec+import Control.Foldl (list, purely)+import Control.Monad.IO.Class (liftIO)+import Data.Vector (fromList)+import Pipes (each, runEffect, (>->))+import qualified Pipes.Prelude as P+import Pipes.Safe (runSafeT)+import SequenceFormats.Eigenstrat (EigenstratIndEntry (..),+ EigenstratSnpEntry (..),+ GenoEntry (..), GenoLine, Sex (..),+ readEigenstrat, writeEigenstrat)+import SequenceFormats.Utils (Chrom (..))+import Test.Hspec spec :: Spec spec = do testReadEigenstrat+ testReadEigenstratCompressed testWriteEigenstrat mockDatEigenstratSnp :: [EigenstratSnpEntry]@@ -44,7 +47,7 @@ fromList [HomRef, Het, Het, HomAlt, HomAlt], fromList [HomAlt, HomAlt, Het, Missing, Het], fromList [HomRef, HomRef, Het, Missing, Missing]]- + testReadEigenstrat :: Spec testReadEigenstrat = describe "readEigenstrat" $ do it "should read the correct eigenstrat file" $ do@@ -52,11 +55,24 @@ esIndFile = "testDat/example.ind" esGenoFile = "testDat/example.eigenstratgeno" (indEntries, esProd) <- runSafeT $ readEigenstrat esGenoFile esSnpFile esIndFile- indEntries `shouldBe` mockDatEigenstratInd + indEntries `shouldBe` mockDatEigenstratInd snpGenoEntries <- runSafeT $ purely P.fold list esProd- (map fst snpGenoEntries) `shouldBe` mockDatEigenstratSnp + (map fst snpGenoEntries) `shouldBe` mockDatEigenstratSnp (map snd snpGenoEntries) `shouldBe` mockDatEigenstratGeno +testReadEigenstratCompressed :: Spec+testReadEigenstratCompressed = describe "readEigenstrat with gzip" $ do+ it "should read the correct eigenstrat file" $ do+ let esSnpFile = "testDat/example.snp.gz"+ esIndFile = "testDat/example.ind"+ esGenoFile = "testDat/example.eigenstratgeno.gz"+ (indEntries, esProd) <- runSafeT $ readEigenstrat esGenoFile esSnpFile esIndFile+ indEntries `shouldBe` mockDatEigenstratInd+ snpGenoEntries <- runSafeT $ purely P.fold list esProd+ (map fst snpGenoEntries) `shouldBe` mockDatEigenstratSnp+ (map snd snpGenoEntries) `shouldBe` mockDatEigenstratGeno++ testWriteEigenstrat :: Spec testWriteEigenstrat = describe "writeEigenstrat" $ do it "should write and read back eigenstrat data correctly" $ do@@ -69,7 +85,7 @@ liftIO . runSafeT . runEffect $ testDatJointProd >-> writeEigenstrat tmpGeno tmpSnp tmpInd mockDatEigenstratInd (indEntries, esProd) <- liftIO . runSafeT $ readEigenstrat tmpGeno tmpSnp tmpInd- indEntries `shouldBe` mockDatEigenstratInd + indEntries `shouldBe` mockDatEigenstratInd snpGenoEntries <- liftIO . runSafeT $ purely P.fold list esProd (map fst snpGenoEntries) `shouldBe` mockDatEigenstratSnp (map snd snpGenoEntries) `shouldBe` mockDatEigenstratGeno
test/SequenceFormats/FreqSumSpec.hs view
@@ -1,15 +1,15 @@ {-# LANGUAGE OverloadedStrings #-} module SequenceFormats.FreqSumSpec (spec) where -import SequenceFormats.FreqSum (readFreqSumFile, printFreqSumFile, FreqSumEntry(..), - FreqSumHeader(..))-import SequenceFormats.Utils (Chrom(..))+import SequenceFormats.FreqSum (FreqSumEntry (..), FreqSumHeader (..),+ printFreqSumFile, readFreqSumFile)+import SequenceFormats.Utils (Chrom (..)) -import Control.Foldl (purely, list)-import Pipes (each, runEffect, (>->))-import qualified Pipes.Prelude as P-import Pipes.Safe (runSafeT)-import Test.Hspec+import Control.Foldl (list, purely)+import Pipes (each, runEffect, (>->))+import qualified Pipes.Prelude as P+import Pipes.Safe (runSafeT)+import Test.Hspec spec :: Spec spec = do@@ -23,9 +23,9 @@ fsEntries_ <- purely P.fold list fsProd_ return (fsHeader_, fsEntries_) it "should read the correct fs header" $- fsHeader `shouldBe` mockDatFsHeader + fsHeader `shouldBe` mockDatFsHeader it "should read the correct fs entries" $- fsEntries `shouldBe` mockDatFsEntries + fsEntries `shouldBe` mockDatFsEntries testPrintFreqSumFile :: Spec testPrintFreqSumFile = describe "printFreqSumFile" $ do@@ -37,11 +37,11 @@ fsEntries_ <- purely P.fold list fsProd_ return (fsHeader_, fsEntries_) it "should read the correct fs header after writing" $- fsHeader `shouldBe` mockDatFsHeader + fsHeader `shouldBe` mockDatFsHeader it "should read the correct fs entries after writing" $- fsEntries `shouldBe` mockDatFsEntries + fsEntries `shouldBe` mockDatFsEntries -mockDatFsHeader :: FreqSumHeader +mockDatFsHeader :: FreqSumHeader mockDatFsHeader = FreqSumHeader names numbers where names = ["SAMPLE0", "SAMPLE1", "SAMPLE2", "SAMPLE3", "SAMPLE4"]
test/SequenceFormats/PileupSpec.hs view
@@ -1,13 +1,14 @@ {-# LANGUAGE OverloadedStrings #-} module SequenceFormats.PileupSpec (spec) where -import SequenceFormats.Pileup (readPileupFromFile, PileupRow(..), Strand(..))-import SequenceFormats.Utils (Chrom(..))+import SequenceFormats.Pileup (PileupRow (..), Strand (..),+ readPileupFromFile)+import SequenceFormats.Utils (Chrom (..)) -import Control.Foldl (purely, list)-import qualified Pipes.Prelude as P-import Pipes.Safe (runSafeT)-import Test.Hspec+import Control.Foldl (list, purely)+import qualified Pipes.Prelude as P+import Pipes.Safe (runSafeT)+import Test.Hspec spec :: Spec spec = testReadPileupFromFile@@ -21,7 +22,7 @@ mockDatPentries :: [PileupRow] mockDatPentries = [ PileupRow (Chrom "1") 1000 'A' ["AAACA", "AAAC", "AAAACCAACA"]- [[f, f, r, f, f], [f, f, r, r], [r, r, f, r, r, r, f, f, r, f]], + [[f, f, r, f, f], [f, f, r, r], [r, r, f, r, r, r, f, f, r, f]], PileupRow (Chrom "1") 2000 'C' ["CCCA", "ACTCC", "CACACCCC"] [[f, f, f, r], [f, f, r, f, f], [f, r, f, r, f, f, f, f]], PileupRow (Chrom "2") 1000 'G' ["GGGGGGGCG", "GGG", "GGGGGG"]
test/SequenceFormats/PlinkSpec.hs view
@@ -4,11 +4,13 @@ import SequenceFormats.Eigenstrat (EigenstratIndEntry (..), EigenstratSnpEntry (..), GenoEntry (..), GenoLine, Sex (..))-import SequenceFormats.Plink (readBimFile, readFamFile,+import SequenceFormats.Plink (PlinkFamEntry (..),+ PlinkPopNameMode (..),+ eigenstratInd2PlinkFam,+ plinkFam2EigenstratInd,+ readBimFile, readFamFile, readPlink, readPlinkBedFile,- writePlink, PlinkFamEntry(..),- plinkFam2EigenstratInd, eigenstratInd2PlinkFam,- PlinkPopNameMode(..))+ writePlink) import SequenceFormats.Utils (Chrom (..)) import Control.Foldl (list, purely)@@ -22,8 +24,10 @@ spec :: Spec spec = do testReadBimFile+ testReadBimFileCompressed testReadFamFile testReadBedFile+ testReadBedFileCompressed testReadPlink testWritePlink testFam2Ind@@ -60,9 +64,15 @@ testReadBimFile :: Spec testReadBimFile = describe "readBimFile" $ it "should read a BIM file correctly" $ do- let esSnpProd = readBimFile "testDat/example.bim"+ let esSnpProd = readBimFile "testDat/example.plink.bim" (runSafeT $ purely P.fold list esSnpProd) `shouldReturn` mockDatEigenstratSnp +testReadBimFileCompressed :: Spec+testReadBimFileCompressed = describe "readBimFile with gzip" $+ it "should read a BIM file correctly" $ do+ let esSnpProd = readBimFile "testDat/example.plink.bim.gz"+ (runSafeT $ purely P.fold list esSnpProd) `shouldReturn` mockDatEigenstratSnp+ testReadFamFile :: Spec testReadFamFile = describe "readFamFile" $ it "should read a FAM file correctly" $ do@@ -77,6 +87,15 @@ purely P.fold list bedProd bedDat `shouldBe` mockDatPlinkBed +testReadBedFileCompressed :: Spec+testReadBedFileCompressed = describe "readBedFile with gzip" $+ it "should read genotypes correctly" $ do+ let fn = "testDat/example.plink.bed.gz"+ bedDat <- runSafeT $ do+ bedProd <- readPlinkBedFile fn 5+ purely P.fold list bedProd+ bedDat `shouldBe` mockDatPlinkBed+ testReadPlink :: Spec testReadPlink = describe "readPlink" $ do it "should read the correct Plink files" $ do@@ -120,7 +139,7 @@ it "should correctly convert with both-option when they are the same" $ do let fam = PlinkFamEntry "Pop1" "SAMPLE0" "0" "0" Female "Pop1" plinkFam2EigenstratInd PlinkPopNameAsBoth fam `shouldBe` EigenstratIndEntry "SAMPLE0" Female "Pop1"- + testInd2Fam :: Spec testInd2Fam = describe "eigenstratInd2PlinkFam" $ do it "should correctly convert with family-option" $ do
test/SequenceFormats/RareAlleleHistogramSpec.hs view
@@ -1,10 +1,12 @@ {-# LANGUAGE OverloadedStrings #-} module SequenceFormats.RareAlleleHistogramSpec (spec) where -import SequenceFormats.RareAlleleHistogram (RareAlleleHistogram(..), readHistogram, writeHistogramFile)+import SequenceFormats.RareAlleleHistogram (RareAlleleHistogram (..),+ readHistogram,+ writeHistogramFile) -import qualified Data.Map as Map-import Test.Hspec+import qualified Data.Map as Map+import Test.Hspec spec :: Spec spec = do
test/SequenceFormats/UtilsSpec.hs view
@@ -1,10 +1,10 @@ {-# LANGUAGE OverloadedStrings #-} module SequenceFormats.UtilsSpec (spec) where -import SequenceFormats.Utils (Chrom(..), SeqFormatException(..))+import SequenceFormats.Utils (Chrom (..), SeqFormatException (..)) -import Control.Exception (evaluate)-import Test.Hspec+import Control.Exception (evaluate)+import Test.Hspec spec :: Spec spec = testChrom@@ -31,5 +31,5 @@ Chrom "X" < Chrom "chrMT" `shouldBe` True specify "chrSSS should throw" $ evaluate (Chrom "chrSSS" < Chrom "chrMT") `shouldThrow` (==SeqFormatException "cannot parse chromosome SSS")- +
test/SequenceFormats/VCFSpec.hs view
@@ -1,18 +1,21 @@ {-# LANGUAGE OverloadedStrings #-} module SequenceFormats.VCFSpec (spec) where -import Control.Foldl (list, purely)-import Pipes.Prelude (fold)-import Pipes.Safe (runSafeT)-import SequenceFormats.FreqSum (FreqSumEntry(..))-import SequenceFormats.Utils (Chrom(..))-import SequenceFormats.VCF (readVCFfromFile, getGenotypes, getDosages,- isTransversionSnp, vcfToFreqSumEntry, isBiallelicSnp, VCFheader(..), VCFentry(..))-import Test.Hspec+import Control.Foldl (list, purely)+import Pipes.Prelude (fold)+import Pipes.Safe (runSafeT)+import SequenceFormats.FreqSum (FreqSumEntry (..))+import SequenceFormats.Utils (Chrom (..))+import SequenceFormats.VCF (VCFentry (..), VCFheader (..),+ getDosages, getGenotypes,+ isBiallelicSnp, isTransversionSnp,+ readVCFfromFile, vcfToFreqSumEntry)+import Test.Hspec spec :: Spec spec = do testReadVCFfromFile+ testReadVCFfromFileCompressed testGetGenotypes testGetDosages testIsTransversionSnp@@ -30,7 +33,23 @@ vcfHc !! 0 `shouldBe` "##fileformat=VCFv4.2" vcfHc !! 18 `shouldBe` "##bcftools_callCommand=call -c -v" it "reads the correct sample names" $- vcfSampleNames vcfH `shouldBe` ["12880A", "12881A", "12883A", "12884A", "12885A"] + vcfSampleNames vcfH `shouldBe` ["12880A", "12881A", "12883A", "12884A", "12885A"]+ it "reads the correct vcf genotype rows" $ do+ vcfRows !! 0 `shouldBe` vcf1+ vcfRows !! 6 `shouldBe` vcf7++testReadVCFfromFileCompressed :: Spec+testReadVCFfromFileCompressed = describe "readVCFfromFile with gzip" $ do+ (vcfH, vcfRows) <- runIO . runSafeT $ do+ (vcfH_, vcfProd_) <- readVCFfromFile "testDat/example.vcf.gz"+ vcfRows_ <- purely fold list vcfProd_+ return (vcfH_, vcfRows_)+ let vcfHc = vcfHeaderComments vcfH+ it "reads the correct header lines" $ do+ vcfHc !! 0 `shouldBe` "##fileformat=VCFv4.2"+ vcfHc !! 18 `shouldBe` "##bcftools_callCommand=call -c -v"+ it "reads the correct sample names" $+ vcfSampleNames vcfH `shouldBe` ["12880A", "12881A", "12883A", "12884A", "12885A"] it "reads the correct vcf genotype rows" $ do vcfRows !! 0 `shouldBe` vcf1 vcfRows !! 6 `shouldBe` vcf7
− testDat/example.bim
@@ -1,7 +0,0 @@-11 rs0000 0.000000 0 A C-11 rs1111 0.001000 100000 A G-11 rs2222 0.002000 200000 A T-11 rs3333 0.003000 300000 C A-11 rs4444 0.004000 400000 G A-11 rs5555 0.005000 500000 T A-11 rs6666 0.006000 600000 G T