diff --git a/sequence-formats.cabal b/sequence-formats.cabal --- a/sequence-formats.cabal +++ b/sequence-formats.cabal @@ -1,78 +1,54 @@ -cabal-version: >=1.10 -name: sequence-formats -version: 1.4.1 -license: GPL-3 -license-file: LICENSE -maintainer: stephan.schiffels@mac.com -author: Stephan Schiffels -synopsis: A package with basic parsing utilities for several Bioinformatic data formats. -description: - Contains utilities to parse and write Eigenstrat, Fasta, FreqSum, VCF and other file formats used in population genetics analyses. -category: Bioinformatics -build-type: Simple -extra-source-files: - README.md - Changelog.md - testDat/example.bim - 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 +name: sequence-formats +version: 1.5.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 +extra-source-files: README.md, + Changelog.md, + testDat/example.bim, + 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 library - exposed-modules: - SequenceFormats.RareAlleleHistogram - SequenceFormats.FreqSum - SequenceFormats.Fasta - SequenceFormats.VCF - SequenceFormats.Eigenstrat - SequenceFormats.Utils - SequenceFormats.Pileup - hs-source-dirs: src - default-language: Haskell2010 - build-depends: - base >=4.7 && <5, - containers >=0.6.0.1, - errors >=2.3.0, - attoparsec >=0.13.2.2, - pipes >=4.3.11, - transformers >=0.5.6.2, - bytestring >=0.10.8.2, - lens-family >=1.2.3, - pipes-bytestring >=2.1.6, - foldl >=1.4.5, - exceptions >=0.10.2, - pipes-safe >=2.3.1, - pipes-attoparsec >=0.5.1.5, - vector >=0.12.0.3 + exposed-modules: SequenceFormats.RareAlleleHistogram, + SequenceFormats.FreqSum, + SequenceFormats.Fasta, + SequenceFormats.VCF, + SequenceFormats.Eigenstrat, + SequenceFormats.Plink, + SequenceFormats.Utils, + SequenceFormats.Pileup + 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 + default-language: Haskell2010 -test-suite sequenceFormatTests - type: exitcode-stdio-1.0 - main-is: Spec.hs - hs-source-dirs: test - other-modules: - SequenceFormats.EigenstratSpec - SequenceFormats.FastaSpec - SequenceFormats.FreqSumSpec - SequenceFormats.RareAlleleHistogramSpec - SequenceFormats.UtilsSpec - SequenceFormats.VCFSpec - SequenceFormats.PileupSpec - default-language: Haskell2010 - build-depends: - base >=4.12.0.0, - sequence-formats -any, - foldl >=1.4.5, - pipes >=4.3.11, - pipes-safe >=2.3.1, - tasty >=1.2.3, - vector >=0.12.0.3, - transformers >=0.5.6.2, - tasty-hunit >=0.10.0.2, - bytestring >=0.10.8.2, - containers >=0.6.0.1, - hspec >=2.7.1 +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.FastaSpec, + SequenceFormats.FreqSumSpec, + SequenceFormats.RareAlleleHistogramSpec, + SequenceFormats.UtilsSpec, + SequenceFormats.PlinkSpec, + SequenceFormats.VCFSpec, + SequenceFormats.PileupSpec + default-language: Haskell2010 diff --git a/src/SequenceFormats/Eigenstrat.hs b/src/SequenceFormats/Eigenstrat.hs --- a/src/SequenceFormats/Eigenstrat.hs +++ b/src/SequenceFormats/Eigenstrat.hs @@ -4,50 +4,63 @@ -} -module SequenceFormats.Eigenstrat (EigenstratSnpEntry(..), EigenstratIndEntry(..), - readEigenstratInd, GenoEntry(..), GenoLine, Sex(..), - readEigenstratSnpStdIn, readEigenstratSnpFile, readBimStdIn, readBimFile, - readEigenstrat, writeEigenstrat, writeEigenstratIndFile, writeEigenstratSnp, writeBim, +module SequenceFormats.Eigenstrat (EigenstratSnpEntry(..), EigenstratIndEntry(..), + readEigenstratInd, GenoEntry(..), GenoLine, Sex(..), + readEigenstratSnpStdIn, readEigenstratSnpFile, + readEigenstrat, writeEigenstrat, writeEigenstratIndFile, writeEigenstratSnp, writeEigenstratGeno) where -import SequenceFormats.Utils (consumeProducer, SeqFormatException(..), - Chrom(..), readFileProd, word) +import SequenceFormats.Utils (Chrom (..), + SeqFormatException (..), + consumeProducer, + readFileProd, word) -import Control.Applicative ((<|>)) -import Control.Exception (throw) -import Control.Monad (void, forM_) -import Control.Monad.Catch (MonadThrow) -import Control.Monad.IO.Class (MonadIO, liftIO) +import Control.Applicative ((<|>)) +import Control.Exception (throw) +import Control.Monad (forM_, void) +import Control.Monad.Catch (MonadThrow) +import Control.Monad.IO.Class (MonadIO, liftIO) import qualified Data.Attoparsec.ByteString.Char8 as A -import Data.Vector (Vector, fromList, toList) -import qualified Data.ByteString.Char8 as B -import Pipes (Producer, Pipe, (>->), for, cat, yield, Consumer) -import Pipes.Safe (MonadSafe) -import qualified Pipes.Safe.Prelude as PS -import qualified Pipes.Prelude as P -import qualified Pipes.ByteString as PB -import System.IO (withFile, IOMode(..), Handle, hPutStrLn) +import qualified Data.ByteString.Char8 as B +import Data.Vector (Vector, fromList, toList) +import Pipes (Consumer, Pipe, Producer, + cat, for, yield, (>->)) +import qualified Pipes.ByteString as PB +import qualified Pipes.Prelude as P +import Pipes.Safe (MonadSafe) +import qualified Pipes.Safe.Prelude as PS +import System.IO (Handle, IOMode (..), + hPutStrLn, withFile) -- |A datatype to represent a single genomic SNP. The constructor arguments are: -- Chromosome, Position, Reference Allele, Alternative Allele. -data EigenstratSnpEntry = EigenstratSnpEntry { - snpChrom :: Chrom, - snpPos :: Int, - snpGeneticPos :: Double, - snpId :: B.ByteString, - snpRef :: Char, - snpAlt :: Char - } deriving (Eq, Show) +data EigenstratSnpEntry = EigenstratSnpEntry + { snpChrom :: Chrom + , snpPos :: Int + , snpGeneticPos :: Double + , snpId :: B.ByteString + , snpRef :: Char + , snpAlt :: Char + } + deriving (Eq, Show) -- |A datatype to represent a single individual. The constructor arguments are: -- Name, Sex and Population Name -data EigenstratIndEntry = EigenstratIndEntry String Sex String deriving (Eq, Show) +data EigenstratIndEntry = EigenstratIndEntry String Sex String + deriving (Eq, Show) -- |A datatype to represent Sex in an Eigenstrat Individual file -data Sex = Male | Female | Unknown deriving (Eq, Show) +data Sex = Male + | Female + | Unknown + deriving (Eq, Show) -- |A datatype to represent the genotype of an individual at a SNP. -data GenoEntry = HomRef | Het | HomAlt | Missing deriving (Eq, Show) +data GenoEntry = HomRef + | Het + | HomAlt + | Missing + deriving (Eq, Show) -- |Vector of the genotypes of all individuals at a single SNP. type GenoLine = Vector GenoEntry @@ -63,17 +76,6 @@ void A.endOfLine return $ EigenstratSnpEntry (Chrom chrom) pos geneticPos snpId_ ref alt -bimParser :: A.Parser EigenstratSnpEntry -bimParser = do - chrom <- word - snpId_ <- A.skipMany1 A.space >> word - geneticPos <- A.skipMany1 A.space >> A.double - pos <- A.skipMany1 A.space >> A.decimal - ref <- A.skipMany1 A.space >> A.satisfy (A.inClass "ACTGN") - alt <- A.skipMany1 A.space >> A.satisfy (A.inClass "ACTGX") - void A.endOfLine - return $ EigenstratSnpEntry (Chrom chrom) pos geneticPos snpId_ ref alt - eigenstratIndParser :: A.Parser EigenstratIndEntry eigenstratIndParser = do A.skipMany A.space @@ -109,7 +111,7 @@ '1' -> return Het '2' -> return HomRef '9' -> return Missing - _ -> error "this should never happen" + _ -> error "this should never happen" where isValidNum c = c == '0' || c == '1' || c == '2' || c == '9' @@ -121,14 +123,6 @@ readEigenstratSnpFile :: (MonadSafe m) => FilePath -> Producer EigenstratSnpEntry m () readEigenstratSnpFile = consumeProducer eigenstratSnpParser . readFileProd --- |Function to read a Bim File from StdIn. Returns a Pipes-Producer over the EigenstratSnpEntries. -readBimStdIn :: (MonadThrow m, MonadIO m) => Producer EigenstratSnpEntry m () -readBimStdIn = consumeProducer bimParser PB.stdin - --- |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 - -- |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 @@ -137,7 +131,7 @@ readEigenstrat genoFile snpFile indFile = do indEntries <- readEigenstratInd indFile let snpProd = readEigenstratSnpFile snpFile - genoProd = consumeProducer eigenstratGenoParser (readFileProd genoFile) >-> + genoProd = consumeProducer eigenstratGenoParser (readFileProd genoFile) >-> validateEigenstratEntries (length indEntries) return (indEntries, P.zip snpProd genoProd) @@ -159,8 +153,8 @@ hPutStrLn h $ name <> "\t" <> sexToStr sex <> "\t" <> popName where sexToStr sex = case sex of - Male -> "M" - Female -> "F" + Male -> "M" + Female -> "F" Unknown -> "U" -- |Function to write an Eigenstrat Snp File. Returns a consumer expecting EigenstratSnpEntries. @@ -174,15 +168,6 @@ in snpLine <> "\n") in toTextPipe >-> snpOutTextConsumer --- |Function to write a Bim file. Returns a consumer expecting EigenstratSnpEntries. -writeBim :: (MonadIO m) => Handle -- ^The Eigenstrat Snp File handle. - -> Consumer EigenstratSnpEntry m () -- ^A consumer to read EigenstratSnpEntries -writeBim snpFileH = - let snpOutTextConsumer = PB.toHandle snpFileH - toTextPipe = P.map (\(EigenstratSnpEntry chrom pos gpos gid ref alt) -> - B.intercalate "\t" [unChrom chrom, gid, B.pack (show gpos), B.pack (show pos), B.singleton ref, B.singleton alt]) - in toTextPipe >-> snpOutTextConsumer - -- |Function to write an Eigentrat Geno File. Returns a consumer expecting Eigenstrat Genolines. writeEigenstratGeno :: (MonadIO m) => Handle -- ^The Genotype file handle -> Consumer GenoLine m () -- ^A consumer to read Genotype entries. @@ -194,9 +179,9 @@ in toTextPipe >-> genoOutTextConsumer where toEigenStratNum c = case c of - HomRef -> 2 :: Int - Het -> 1 - HomAlt -> 0 + HomRef -> 2 :: Int + Het -> 1 + HomAlt -> 0 Missing -> 9 -- |Function to write an Eigenstrat Database. Returns a consumer expecting joint Snp- and Genotype lines. diff --git a/src/SequenceFormats/Plink.hs b/src/SequenceFormats/Plink.hs new file mode 100644 --- /dev/null +++ b/src/SequenceFormats/Plink.hs @@ -0,0 +1,127 @@ +{-# LANGUAGE BinaryLiterals #-} +{-# LANGUAGE OverloadedStrings #-} +module SequenceFormats.Plink (readBimStdIn, readBimFile, writeBim, readFamFile, readPlinkBedFile) where + +import SequenceFormats.Eigenstrat (EigenstratIndEntry (..), + EigenstratSnpEntry (..), + GenoEntry (..), GenoLine, + Sex (..)) +import SequenceFormats.Utils (Chrom (..), consumeProducer, + readFileProd, word) + +import Control.Applicative ((<|>)) +import Control.Monad (void) +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 as AB +import qualified Data.Attoparsec.ByteString.Char8 as A +import Data.Bits (shiftR, (.&.)) +import qualified Data.ByteString as BB +import qualified Data.ByteString.Char8 as B +import Data.Vector (fromList) +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, runSafeT) +import qualified Pipes.Safe.Prelude +import System.IO (Handle, IOMode (..), + withFile) + + +bimParser :: A.Parser EigenstratSnpEntry +bimParser = do + chrom <- word + snpId_ <- A.skipMany1 A.space >> word + geneticPos <- A.skipMany1 A.space >> A.double + pos <- A.skipMany1 A.space >> A.decimal + ref <- A.skipMany1 A.space >> A.satisfy (A.inClass "ACTGN") + alt <- A.skipMany1 A.space >> A.satisfy (A.inClass "ACTGX") + void A.endOfLine + return $ EigenstratSnpEntry (Chrom chrom) pos geneticPos snpId_ ref alt + +famParser :: A.Parser EigenstratIndEntry +famParser = do + A.skipMany A.space + pop <- word + ind <- A.skipMany1 A.space >> word + _ <- A.skipMany1 A.space >> A.decimal + _ <- A.skipMany1 A.space >> A.decimal + sex <- A.skipMany1 A.space >> parseSex + _ <- A.skipMany1 A.space >> word + void A.endOfLine + return $ EigenstratIndEntry (B.unpack ind) sex (B.unpack pop) + where + parseSex = parseMale <|> parseFemale <|> parseUnknown + parseMale = A.char '1' >> return Male + parseFemale = A.char '2' >> return Female + parseUnknown = A.anyChar >> return Unknown + +bedHeaderParser :: AB.Parser () +bedHeaderParser = do + void $ AB.word8 0b01101100 -- magic number I for BED files + void $ AB.word8 0b00011011 -- magic number II for BED files + void $ AB.word8 0b00000001 -- we can only parse SNP-major order + +bedGenotypeParser :: Int -> AB.Parser GenoLine +bedGenotypeParser nrInds = do + let nrBytes = if nrInds `rem` 4 == 0 then nrInds `quot` 4 else (nrInds `quot` 4) + 1 + bytes <- BB.unpack <$> AB.take nrBytes + let indBitPairs = concatMap getBitPairs bytes + return . fromList . take nrInds . map bitPairToGenotype $ indBitPairs + where + getBitPairs byte = map (0b00000011 .&.) [byte, shiftR byte 2, shiftR byte 4, shiftR byte 6] + bitPairToGenotype 0b00000000 = HomRef + bitPairToGenotype 0b00000010 = Het + bitPairToGenotype 0b00000011 = HomAlt + bitPairToGenotype 0b00000001 = Missing + bitPairToGenotype _ = error "This should never happen" + +readPlinkBedProd :: (MonadThrow m) => Int -> Producer B.ByteString m () -> m (Producer GenoLine m ()) +readPlinkBedProd nrInds prod = do + (res, rest) <- runStateT (parse bedHeaderParser) prod + _ <- case res of + Nothing -> throwM $ ParsingError [] "Bed file exhausted prematurely" + Just (Left e) -> throwM e + Just (Right h) -> return h + return $ consumeProducer (bedGenotypeParser nrInds) rest + +-- |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 (Pipes.Safe.Prelude.withFile file ReadMode PB.fromHandle) + +-- |Function to read a Bim File from StdIn. Returns a Pipes-Producer over the EigenstratSnpEntries. +readBimStdIn :: (MonadThrow m, MonadIO m) => Producer EigenstratSnpEntry m () +readBimStdIn = consumeProducer bimParser PB.stdin + +-- |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 + +-- |Function to read a Plink fam file. Returns the Eigenstrat Individual Entries as list. +readFamFile :: (MonadIO m) => FilePath -> m [EigenstratIndEntry] +readFamFile fn = + liftIO . withFile fn ReadMode $ \handle -> + P.toListM $ consumeProducer famParser (PB.fromHandle handle) + +-- |Function to read a full Plink dataset from files. Returns a pair of the Plink Individual Entries, and a joint Producer over the snp entries and the genotypes. +readPlink :: (MonadSafe m) => FilePath -- ^The Bed file + -> FilePath -- ^The Bim File + -> FilePath -- ^The Fam file + -> m ([EigenstratIndEntry], Producer (EigenstratSnpEntry, GenoLine) m ()) -- The return pair of individual entries and a joint Snp/Geno Producer. +readPlink bedFile bimFile famFile = do + indEntries <- readFamFile famFile + let nrInds = length indEntries + snpProd = readBimFile bimFile + genoProd <- readPlinkBedFile bedFile nrInds + return (indEntries, P.zip snpProd genoProd) + +-- |Function to write a Bim file. Returns a consumer expecting EigenstratSnpEntries. +writeBim :: (MonadIO m) => Handle -- ^The Eigenstrat Snp File handle. + -> Consumer EigenstratSnpEntry m () -- ^A consumer to read EigenstratSnpEntries +writeBim snpFileH = + let snpOutTextConsumer = PB.toHandle snpFileH + toTextPipe = P.map (\(EigenstratSnpEntry chrom pos gpos gid ref alt) -> + B.intercalate "\t" [unChrom chrom, gid, B.pack (show gpos), B.pack (show pos), B.singleton ref, B.singleton alt]) + in toTextPipe >-> snpOutTextConsumer diff --git a/test/SequenceFormats/EigenstratSpec.hs b/test/SequenceFormats/EigenstratSpec.hs --- a/test/SequenceFormats/EigenstratSpec.hs +++ b/test/SequenceFormats/EigenstratSpec.hs @@ -7,14 +7,13 @@ import Pipes (each, runEffect, (>->)) import qualified Pipes.Prelude as P import Pipes.Safe (runSafeT) -import SequenceFormats.Eigenstrat (readEigenstrat, writeEigenstrat, readBimFile, +import SequenceFormats.Eigenstrat (readEigenstrat, writeEigenstrat, EigenstratSnpEntry(..), EigenstratIndEntry(..), GenoLine, Sex(..), GenoEntry(..)) import SequenceFormats.Utils (Chrom(..)) import Test.Hspec spec :: Spec spec = do - testReadBimFile testReadEigenstrat testWriteEigenstrat @@ -27,12 +26,6 @@ EigenstratSnpEntry (Chrom "11") 400000 0.004000 "rs4444" 'G' 'A', EigenstratSnpEntry (Chrom "11") 500000 0.005000 "rs5555" 'T' 'A', EigenstratSnpEntry (Chrom "11") 600000 0.006000 "rs6666" 'G' 'T'] - -testReadBimFile :: Spec -testReadBimFile = describe "readBimFile" $ - it "should read a BIM file correctly" $ do - let esSnpProd = readBimFile "testDat/example.bim" - (runSafeT $ purely P.fold list esSnpProd) `shouldReturn` mockDatEigenstratSnp mockDatEigenstratInd :: [EigenstratIndEntry] mockDatEigenstratInd = [ diff --git a/test/SequenceFormats/PlinkSpec.hs b/test/SequenceFormats/PlinkSpec.hs new file mode 100644 --- /dev/null +++ b/test/SequenceFormats/PlinkSpec.hs @@ -0,0 +1,67 @@ +{-# LANGUAGE OverloadedStrings #-} +module SequenceFormats.PlinkSpec (spec) where + +import SequenceFormats.Eigenstrat (EigenstratIndEntry (..), + EigenstratSnpEntry (..), Sex (..), GenoLine, GenoEntry(..)) +import SequenceFormats.Plink (readBimFile, readFamFile, readPlinkBedFile) +import SequenceFormats.Utils (Chrom (..)) + +import Control.Foldl (list, purely) +import Data.Vector (fromList) +import qualified Pipes.Prelude as P +import Pipes.Safe (runSafeT) +import Test.Hspec + +spec :: Spec +spec = do + testReadBimFile + testReadFamFile + testReadBedFile + +mockDatEigenstratSnp :: [EigenstratSnpEntry] +mockDatEigenstratSnp = [ + EigenstratSnpEntry (Chrom "11") 0 0.000000 "rs0000" 'A' 'C', + EigenstratSnpEntry (Chrom "11") 100000 0.001000 "rs1111" 'A' 'G', + EigenstratSnpEntry (Chrom "11") 200000 0.002000 "rs2222" 'A' 'T', + EigenstratSnpEntry (Chrom "11") 300000 0.003000 "rs3333" 'C' 'A', + EigenstratSnpEntry (Chrom "11") 400000 0.004000 "rs4444" 'G' 'A', + EigenstratSnpEntry (Chrom "11") 500000 0.005000 "rs5555" 'T' 'A', + EigenstratSnpEntry (Chrom "11") 600000 0.006000 "rs6666" 'G' 'T'] + +mockDatEigenstratInd :: [EigenstratIndEntry] +mockDatEigenstratInd = [ + EigenstratIndEntry "SAMPLE0" Female "1", + EigenstratIndEntry "SAMPLE1" Male "2", + EigenstratIndEntry "SAMPLE2" Female "3", + EigenstratIndEntry "SAMPLE3" Male "4", + EigenstratIndEntry "SAMPLE4" Female "5"] + +mockDatPlinkBed :: [GenoLine] +mockDatPlinkBed = [ + fromList [Het, Het, Het, HomAlt, HomAlt], + fromList [HomAlt, Het, HomRef, Het, HomRef], + fromList [HomRef, Het, Het, HomAlt, Het], + fromList [HomAlt, HomAlt, Het, HomRef, HomRef], + fromList [HomRef, Het, Het, HomAlt, HomAlt], + fromList [HomAlt, HomAlt, Het, Het, Het], + fromList [HomRef, HomRef, Het, Het, HomAlt]] + +testReadBimFile :: Spec +testReadBimFile = describe "readBimFile" $ + it "should read a BIM file correctly" $ do + let esSnpProd = readBimFile "testDat/example.bim" + (runSafeT $ purely P.fold list esSnpProd) `shouldReturn` mockDatEigenstratSnp + +testReadFamFile :: Spec +testReadFamFile = describe "readFamFile" $ + it "should read a FAM file correctly" $ do + readFamFile "testDat/example.fam" `shouldReturn` mockDatEigenstratInd + +testReadBedFile :: Spec +testReadBedFile = describe "readBedFile" $ + it "should read genotypes correctly" $ do + let fn = "testDat/example.plink.bed" + bedDat <- runSafeT $ do + bedProd <- readPlinkBedFile fn 5 + purely P.fold list bedProd + bedDat `shouldBe` mockDatPlinkBed