diff --git a/Changelog.md b/Changelog.md --- a/Changelog.md +++ b/Changelog.md @@ -1,5 +1,8 @@ # Changelog - +- V 1.12.1.0: Clearer errors for broken input files. Reading a Plink bed file now throws a `SeqFormatException` naming the file if it is too short for the header, has wrong magic bytes (with a hint if it looks gzip-compressed), is not in SNP-major mode, or ends with an incomplete SNP record (e.g. a fam file with the wrong number of individuals). Previously this gave cryptic attoparsec errors such as `Failed reading: satisfy`. Reading gzip files now throws a `SeqFormatException` for invalid gzip data (previously `ZlibException (-3)`), and for truncated gzip files, which were previously accepted without error, possibly silently losing data. +- V 1.12.0.0: Fixed reading of multi-member gzip files, such as BGZF-compressed VCFs written by bgzip, bcftools or GATK. Previously only the first gzip member was read, leading to parsing errors or silently truncated data. Breaking change: `vcfToFreqSumEntry` now returns `Maybe FreqSumEntry`, with `Nothing` for indels, multi-allelic sites and non-nucleotide alleles (instead of throwing an exception at the first indel). Writing Eigenstrat SNP and Plink BIM files now throws an exception on empty SNP IDs or IDs containing whitespace, as these cannot be read back. Clients converting from VCF or FreqSum need to provide IDs for sites without one (e.g. `Chrom_Pos`). +- V 1.11.0.4: improved parsing error messages. +- V 1.11.0.3: added dots to allowed characters in Plink files - V 1.11.0.2: exposed parseSex from Eigenstrat - V 1.11.0.1: Allowing missing alternative alleles when converting VCF to FreqSum, and improved error messaging. - V 1.11.0.0: Added support for writing of VCF files, including gzipping. Made some breaking API changes on top, for example @@ -50,23 +53,3 @@ - V 1.1.5: Fixed VCF parser: Now breaks if lines end prematurely - V 1.1.4.2: Exporting readVCFfromProd - V 1.1.4.1: First entry in the Changelog. Added Haddock documentation to all modules and prepare for releasing on Hackage. - - - - - - - - - - - - - - - - - - - - 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.11.0.2 +version: 1.12.1.0 license: GPL-3 license-file: LICENSE maintainer: stephan.schiffels@mac.com @@ -32,10 +32,16 @@ testDat/example.plink.bim testDat/example.plink.bim.gz testDat/example.plink.fam + testDat/example.short.plink.bed + testDat/example.wrongmagic.plink.bed + testDat/example.gzipped.plink.bed + testDat/example.truncated.plink.bed + testDat/example.truncated.plink.bed.gz testDat/example.snp testDat/example.snp.gz testDat/example.vcf testDat/example.vcf.gz + testDat/example.multimember.vcf.gz library exposed-modules: @@ -67,7 +73,6 @@ pipes-safe >=2.3.5, pipes-attoparsec >=0.6.0, vector >=0.13.1.0, - pipes-zlib >=0.4.4.2, streaming-commons >=0.2.2.6 test-suite sequenceFormatTests diff --git a/src/SequenceFormats/Eigenstrat.hs b/src/SequenceFormats/Eigenstrat.hs --- a/src/SequenceFormats/Eigenstrat.hs +++ b/src/SequenceFormats/Eigenstrat.hs @@ -20,12 +20,13 @@ import Control.Applicative ((<|>)) import Control.Exception (throw) -import Control.Monad (forM_, void) -import Control.Monad.Catch (MonadThrow) +import Control.Monad (forM_, void, when) +import Control.Monad.Catch (MonadThrow, throwM) import Control.Monad.IO.Class (MonadIO, liftIO) 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 Data.List (isSuffixOf) import qualified Data.Streaming.Zlib as Z import Data.Vector (Vector, fromList, toList) @@ -177,10 +178,13 @@ return $ gzipConsumer def snpFileH else return $ PB.toHandle snpFileH - let toTextPipe = P.map (\(EigenstratSnpEntry chrom pos gpos gid ref alt) -> + let toTextPipe = P.mapM (\(EigenstratSnpEntry chrom pos gpos gid ref alt) -> do + when (B.null gid || B.any isSpace gid) . throwM . SeqFormatException $ + "invalid SNP ID " ++ show gid ++ " at " ++ show chrom ++ ":" ++ show pos ++ + ". SNP IDs must be non-empty and must not contain whitespace" let snpLine = B.intercalate "\t" [gid, unChrom chrom, B.pack (show gpos), B.pack (show pos), B.singleton ref, B.singleton alt] - in snpLine <> "\n") + return $ snpLine <> "\n") toTextPipe >-> snpOutTextConsumer -- |Function to write an Eigentrat Geno File. Returns a consumer expecting Eigenstrat Genolines. diff --git a/src/SequenceFormats/Plink.hs b/src/SequenceFormats/Plink.hs --- a/src/SequenceFormats/Plink.hs +++ b/src/SequenceFormats/Plink.hs @@ -16,35 +16,39 @@ EigenstratSnpEntry (..), GenoEntry (..), GenoLine, Sex (..)) -import SequenceFormats.Utils (Chrom (..), consumeProducer, +import SequenceFormats.Utils (Chrom (..), + SeqFormatException (..), + consumeProducer, deflateFinaliser, gzipConsumer, readFileProdCheckCompress, word, writeFromPopper) import Control.Applicative ((<|>)) -import Control.Monad (forM_, void) +import Control.Monad (forM_, void, when) 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 import Data.Bits (shiftL, shiftR, (.&.), (.|.)) import qualified Data.ByteString as BB import qualified Data.ByteString.Char8 as B +import Data.Char (isSpace) import Data.List (isSuffixOf) import qualified Data.Streaming.Zlib as Z import Data.Vector (fromList, toList) import Data.Word (Word8) import Pipes (Consumer, Producer, (>->)) -import Pipes.Attoparsec (ParsingError (..), parse) +import Lens.Family2 (view) +import Pipes.Attoparsec (parsed) import qualified Pipes.ByteString as PB import qualified Pipes.Prelude as P import Pipes.Safe (MonadSafe, register) import qualified Pipes.Safe.Prelude as PS import System.IO (IOMode (..), withFile) +import Text.Printf (printf) -- see https://www.cog-genomics.org/plink/2.0/formats#fam data PlinkFamEntry = PlinkFamEntry { @@ -64,19 +68,21 @@ 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 "ACTGNX01234") - alt <- A.skipMany1 A.space >> A.satisfy (A.inClass "ACTGNX01234") + ref <- A.skipMany1 A.space >> A.satisfy (A.inClass "ACTGNX01234.") + alt <- A.skipMany1 A.space >> A.satisfy (A.inClass "ACTGNX01234.") void A.endOfLine - let refConvert = convertNum ref - altConvert = convertNum alt + let refConvert = convertChar ref + altConvert = convertChar alt return $ EigenstratSnpEntry (Chrom chrom) pos geneticPos snpId_ refConvert altConvert where - convertNum '0' = 'N' - convertNum '1' = 'A' - convertNum '2' = 'C' - convertNum '3' = 'G' - convertNum '4' = 'T' - convertNum x = x + convertChar '0' = 'N' + convertChar '1' = 'A' + convertChar '2' = 'C' + convertChar '3' = 'G' + convertChar '4' = 'T' + convertChar 'X' = 'N' + convertChar '.' = 'N' + convertChar x = x famParser :: A.Parser PlinkFamEntry famParser = do @@ -111,16 +117,33 @@ PlinkPopNameAsPhenotype -> PlinkFamEntry "DummyFamily" indId "0" "0" sex popName PlinkPopNameAsBoth -> PlinkFamEntry popName indId "0" "0" sex popName -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 +-- |Checks the three magic bytes at the start of a bed file, see +-- https://www.cog-genomics.org/plink/1.9/formats#bed +checkBedHeader :: (MonadThrow m) => FilePath -> B.ByteString -> m () +checkBedHeader file header + | BB.length header < 3 = throwM . SeqFormatException $ + "Plink bed file " ++ file ++ " is too short (" ++ show (BB.length header) ++ + " bytes) to contain the 3-byte bed header. The file seems to be empty or truncated" + | BB.take 2 header /= BB.pack [0b01101100, 0b00011011] = throwM . SeqFormatException $ + "Plink bed file " ++ file ++ " does not start with the magic bytes 0x6c 0x1b (found " ++ + showBytes (BB.take 2 header) ++ "). It does not seem to be a Plink bed file" ++ gzipHint + | BB.index header 2 /= 0b00000001 = throwM . SeqFormatException $ + "Plink bed file " ++ file ++ " is not in SNP-major mode (third byte is " ++ + showBytes (BB.drop 2 header) ++ " instead of 0x01). Only SNP-major bed files are supported" + | otherwise = return () + where + showBytes = unwords . map (\b -> printf "0x%02x" b) . BB.unpack + gzipHint = if BB.take 2 header == BB.pack [0x1f, 0x8b] && not (".gz" `isSuffixOf` file) + then ". It looks gzip-compressed, but its name does not end in .gz" + else "" +-- |The number of bytes per SNP in a bed file, with 2 bits per individual. +bedRecordBytes :: Int -> Int +bedRecordBytes nrInds = (nrInds + 3) `quot` 4 + 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 + bytes <- BB.unpack <$> AB.take (bedRecordBytes nrInds) let indBitPairs = concatMap getBitPairs bytes return . fromList . take nrInds . map bitPairToGenotype $ indBitPairs where @@ -131,18 +154,23 @@ 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 +readPlinkBedProd :: (MonadThrow m) => FilePath -> Int -> Producer B.ByteString m () -> m (Producer GenoLine m ()) +readPlinkBedProd file nrInds prod = do + (headerChunks, rest) <- P.toListM' $ view (PB.splitAt (3 :: Int)) prod + checkBedHeader file (BB.concat headerChunks) + return $ parsed (bedGenotypeParser nrInds) rest >>= either reportTruncated return + where + -- the genotype parser can only fail if the input ends in the middle of a SNP record + reportTruncated (_, leftovers) = do + nrLeftoverBytes <- lift $ P.fold (\n c -> n + BB.length c) 0 id leftovers + throwM . SeqFormatException $ "Plink bed file " ++ file ++ " ends with an incomplete SNP " ++ + "record (" ++ show nrLeftoverBytes ++ " bytes left over, but each SNP takes " ++ + show (bedRecordBytes nrInds) ++ " bytes for " ++ show nrInds ++ " individuals). " ++ + "Either the bed file is truncated or it does not match the number of individuals in the fam file" -- |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 . readFileProdCheckCompress $ file +readPlinkBedFile file nrInds = readPlinkBedProd file 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 () @@ -181,10 +209,13 @@ return $ gzipConsumer def bimFileH else return $ PB.toHandle bimFileH - let toTextPipe = P.map (\(EigenstratSnpEntry chrom pos gpos gid ref alt) -> + let toTextPipe = P.mapM (\(EigenstratSnpEntry chrom pos gpos gid ref alt) -> do + when (B.null gid || B.any isSpace gid) . throwM . SeqFormatException $ + "invalid SNP ID " ++ show gid ++ " at " ++ show chrom ++ ":" ++ show pos ++ + ". SNP IDs must be non-empty and must not contain whitespace" let bimLine = B.intercalate "\t" [unChrom chrom, gid, B.pack (show gpos), B.pack (show pos), B.singleton ref, B.singleton alt] - in bimLine <> "\n") + return $ bimLine <> "\n") toTextPipe >-> bimOutTextConsumer -- |Function to write a Plink Fam file. diff --git a/src/SequenceFormats/Utils.hs b/src/SequenceFormats/Utils.hs --- a/src/SequenceFormats/Utils.hs +++ b/src/SequenceFormats/Utils.hs @@ -4,11 +4,13 @@ module SequenceFormats.Utils (liftParsingErrors, consumeProducer, readFileProd, readFileProdCheckCompress, + decompressMultiMember, SeqFormatException(..), deflateFinaliser, Chrom(..), word, gzipConsumer, writeFromPopper, Z.Deflate) where import Control.Error (readErr) import Control.Exception (Exception, throw, throwIO) +import Control.Monad (unless) import Control.Monad.Catch (MonadThrow, throwM) import Control.Monad.IO.Class (MonadIO, liftIO) import Control.Monad.Trans.Class (lift) @@ -18,18 +20,20 @@ import Data.List (isSuffixOf) import qualified Data.Streaming.Zlib as Z import Pipes (Consumer, Producer, await, - next) + next, yield) 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 (Handle, IOMode (..)) -- |An exception type for parsing BioInformatic file formats. data SeqFormatException = SeqFormatException String - deriving (Show, Eq) + deriving (Eq) +instance Show SeqFormatException where + show (SeqFormatException msg) = "SeqFormatException: " ++ msg + instance Exception SeqFormatException -- |A wrapper datatype for Chromosome names. @@ -68,7 +72,8 @@ x <- lift $ next restProd case x of Right (chunk, _) -> do - let msg' = "Error while parsing: " <> msg <> ". Error occurred when trying to parse this chunk: " ++ show chunk + let firstLine = B.unpack . B.takeWhile (/= '\n') $ chunk + msg' = "Error while parsing: " <> msg <> ". Offending line: " ++ firstLine throwM $ SeqFormatException msg' Left _ -> error "should not happen" Right () -> return () @@ -82,8 +87,55 @@ readFileProdCheckCompress :: (PS.MonadSafe m) => FilePath -> Producer B.ByteString m () readFileProdCheckCompress f = - let decompressFunc = if ".gz" `isSuffixOf` f then decompress else id + let decompressFunc = if ".gz" `isSuffixOf` f then decompressGzip ("gzip file " ++ f) else id in decompressFunc $ PS.withFile f ReadMode PB.fromHandle + +-- |Decompresses a gzip stream that may consist of multiple concatenated gzip members, as is the +-- case for BGZF files written by bgzip, bcftools or GATK. Pipes.GZip.decompress on its own stops +-- after the first member and silently drops the rest of the input. Throws a SeqFormatException +-- if the input is not valid gzip data, or if it ends in the middle of a gzip member (e.g. a +-- truncated file), which Pipes.GZip.decompress would silently accept. +decompressMultiMember :: (MonadIO m) => Producer B.ByteString m r -> Producer B.ByteString m r +decompressMultiMember = decompressGzip "gzip stream" + +-- |Like decompressMultiMember, but takes a description of the input (e.g. the file name) for +-- error messages. +decompressGzip :: (MonadIO m) => String -> Producer B.ByteString m r -> Producer B.ByteString m r +decompressGzip descr = newMember True + where + newMember isFirst prod = liftIO (Z.initInflate (Z.WindowBits 31)) >>= go isFirst False prod + -- hasInput tracks whether the current member has received any bytes yet, so that we can + -- tell a clean end of input (after a completed member) from a truncated member. + go isFirst hasInput prod inf = do + res <- lift (next prod) + case res of + Left r -> do + complete <- liftIO (Z.isCompleteInflate inf) + if complete || (not hasInput && not isFirst) then return r else + liftIO . throwIO . SeqFormatException $ if hasInput + then descr ++ " ended unexpectedly. The file seems to be truncated" + else descr ++ " is empty, which is not valid gzip" + Right (bs, prod') + | B.null bs -> go isFirst hasInput prod' inf + | otherwise -> do + popper <- liftIO (Z.feedInflate inf bs) + yieldPopper popper + rest <- liftIO (Z.flushInflate inf) + unless (B.null rest) (yield rest) + complete <- liftIO (Z.isCompleteInflate inf) + if complete then do + leftover <- liftIO (Z.getUnusedInflate inf) + newMember False (yield leftover >> prod') + else + go isFirst True prod' inf + yieldPopper popper = do + popRes <- liftIO popper + case popRes of + Z.PRDone -> return () + Z.PRNext bs -> yield bs >> yieldPopper popper + Z.PRError (Z.ZlibException code) -> liftIO . throwIO . SeqFormatException $ + "could not decompress " ++ descr ++ " (zlib error code " ++ show code ++ + "). The file seems to be corrupt or not gzip-compressed" word :: A.Parser B.ByteString word = A.takeTill isSpace diff --git a/src/SequenceFormats/VCF.hs b/src/SequenceFormats/VCF.hs --- a/src/SequenceFormats/VCF.hs +++ b/src/SequenceFormats/VCF.hs @@ -28,13 +28,14 @@ import Control.Applicative ((<|>)) import Control.Error (atErr) -import Control.Monad (forM, unless, void) +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.Char8 as A import qualified Data.ByteString.Char8 as B +import Data.Char (isAlpha) import Data.List (isSuffixOf) import Data.Maybe (fromMaybe) import qualified Data.Streaming.Zlib as Z @@ -118,7 +119,7 @@ parseFilter = (parseDot *> pure Nothing) <|> (Just <$> word) parseInfoFields = (parseDot *> pure []) <|> (parseInfoField `A.sepBy1` A.char ';') parseInfoField = A.takeTill (\c -> c == ';' || c == '\t') - parseFormatStringsAndGenotypes = (\f g -> Just (f, g)) <$> parseFormatStrings <* sp <*> parseGenotypeInfos + parseFormatStringsAndGenotypes = (\f g -> Just (f, g)) <$> parseFormatStrings <* sp <*> parseGenotypeInfos parseFormatStrings = parseFormatString `A.sepBy1` A.char ':' parseFormatString = A.takeTill (\c -> c == ':' || c == '\t') parseGenotypeInfos = parseGenotype `A.sepBy1` sp @@ -173,21 +174,18 @@ ["1", "1"] -> return $ Just (2, 2) _ -> return Nothing --- |Converts a VCFentry to the simpler FreqSum format -vcfToFreqSumEntry :: (MonadThrow m) => VCFentry -> m FreqSumEntry -vcfToFreqSumEntry vcfEntry = do - unless (B.length (vcfRef vcfEntry) == 1) . throwM . SeqFormatException $ - "multi-site reference allele at " ++ show vcfEntry - alt <- case vcfAlt vcfEntry of - [] -> return 'N' - (a:_) -> if B.length a /= 1 - then - throwM . SeqFormatException $ "multi-site alternative allele at " ++ show vcfEntry - else - return $ B.head a - let ref = B.head (vcfRef vcfEntry) - dosages <- getDosages vcfEntry - return $ FreqSumEntry (vcfChrom vcfEntry) (vcfPos vcfEntry) (vcfId vcfEntry) Nothing ref alt dosages +-- |Converts a VCFentry to the simpler FreqSum format. Returns Nothing for sites that cannot be represented +-- as a biallelic SNP: indels and other multi-base alleles, multi-allelic sites and non-nucleotide alleles such +-- as the spanning deletion allele *. Sites without an alternative allele are kept, with alternative allele N. +vcfToFreqSumEntry :: (MonadThrow m) => VCFentry -> m (Maybe FreqSumEntry) +vcfToFreqSumEntry vcfEntry = case (vcfRef vcfEntry, vcfAlt vcfEntry) of + (ref, []) | isSingleBase ref -> Just <$> makeEntry (B.head ref) 'N' + (ref, [alt]) | isSingleBase ref && isSingleBase alt -> Just <$> makeEntry (B.head ref) (B.head alt) + _ -> return Nothing + where + isSingleBase a = B.length a == 1 && isAlpha (B.head a) + makeEntry ref alt = FreqSumEntry (vcfChrom vcfEntry) (vcfPos vcfEntry) (vcfId vcfEntry) Nothing ref alt <$> + getDosages vcfEntry printVCFtoStdOut :: (MonadIO m) => VCFheader -> Consumer VCFentry m () printVCFtoStdOut vcfh = do diff --git a/test/SequenceFormats/EigenstratSpec.hs b/test/SequenceFormats/EigenstratSpec.hs --- a/test/SequenceFormats/EigenstratSpec.hs +++ b/test/SequenceFormats/EigenstratSpec.hs @@ -3,6 +3,7 @@ import Control.Foldl (list, purely) import Control.Monad.IO.Class (liftIO) +import Data.List (isPrefixOf) import Data.Vector (fromList) import Pipes (each, runEffect, (>->)) import qualified Pipes.Prelude as P @@ -10,8 +11,10 @@ import SequenceFormats.Eigenstrat (EigenstratIndEntry (..), EigenstratSnpEntry (..), GenoEntry (..), GenoLine, Sex (..), - readEigenstrat, writeEigenstrat) -import SequenceFormats.Utils (Chrom (..)) + readEigenstrat, writeEigenstrat, + writeEigenstratSnp) +import SequenceFormats.Utils (Chrom (..), + SeqFormatException (..)) import Test.Hspec spec :: Spec @@ -20,6 +23,7 @@ testReadEigenstratCompressed testWriteEigenstrat testWriteEigenstratCompressed + testWriteEigenstratSnpInvalidId mockDatEigenstratSnp :: [EigenstratSnpEntry] mockDatEigenstratSnp = [ @@ -107,3 +111,10 @@ snpGenoEntries <- liftIO . runSafeT $ purely P.fold list esProd (map fst snpGenoEntries) `shouldBe` mockDatEigenstratSnp (map snd snpGenoEntries) `shouldBe` mockDatEigenstratGeno + +testWriteEigenstratSnpInvalidId :: Spec +testWriteEigenstratSnpInvalidId = describe "writeEigenstratSnp" $ + it "should throw on empty SNP IDs" $ do + let badSnps = [EigenstratSnpEntry (Chrom "11") 0 0.0 "rs0000" 'A' 'C', EigenstratSnpEntry (Chrom "11") 100000 0.001 "" 'A' 'G'] + runSafeT (runEffect (each badSnps >-> writeEigenstratSnp "/tmp/invalidIdTest.snp")) `shouldThrow` + (\(SeqFormatException msg) -> "invalid SNP ID \"\" at 11:100000" `isPrefixOf` msg) diff --git a/test/SequenceFormats/PlinkSpec.hs b/test/SequenceFormats/PlinkSpec.hs --- a/test/SequenceFormats/PlinkSpec.hs +++ b/test/SequenceFormats/PlinkSpec.hs @@ -10,11 +10,13 @@ plinkFam2EigenstratInd, readBimFile, readFamFile, readPlink, readPlinkBedFile, - writePlink) -import SequenceFormats.Utils (Chrom (..)) + writeBim, writePlink) +import SequenceFormats.Utils (Chrom (..), + SeqFormatException (..)) import Control.Foldl (list, purely) import Control.Monad.IO.Class (liftIO) +import Data.List (isInfixOf, isPrefixOf) import Data.Vector (fromList) import Pipes (each, runEffect, (>->)) import qualified Pipes.Prelude as P @@ -28,8 +30,10 @@ testReadFamFile testReadBedFile testReadBedFileCompressed + testReadBedFileInvalid testReadPlink testWritePlink + testWriteBimInvalidId testWritePlinkCompressed testFam2Ind testInd2Fam @@ -97,6 +101,26 @@ purely P.fold list bedProd bedDat `shouldBe` mockDatPlinkBed +testReadBedFileInvalid :: Spec +testReadBedFileInvalid = describe "readBedFile with invalid files" $ do + let readBed fn = runSafeT $ readPlinkBedFile fn 5 >>= P.length + throwsWith substr (SeqFormatException msg) = substr `isInfixOf` msg + it "should report a file too short for the header" $ + readBed "testDat/example.short.plink.bed" `shouldThrow` + throwsWith "is too short (2 bytes) to contain the 3-byte bed header" + it "should report wrong magic bytes" $ + readBed "testDat/example.wrongmagic.plink.bed" `shouldThrow` + throwsWith "does not start with the magic bytes 0x6c 0x1b (found 0x74 0x68)" + it "should report a gzipped file without .gz ending" $ + readBed "testDat/example.gzipped.plink.bed" `shouldThrow` + throwsWith "It looks gzip-compressed, but its name does not end in .gz" + it "should report an incomplete SNP record" $ + readBed "testDat/example.truncated.plink.bed" `shouldThrow` + throwsWith "ends with an incomplete SNP record (1 bytes left over, but each SNP takes 2 bytes for 5 individuals)" + it "should report a truncated gzip file" $ + readBed "testDat/example.truncated.plink.bed.gz" `shouldThrow` + throwsWith "gzip file testDat/example.truncated.plink.bed.gz ended unexpectedly" + testReadPlink :: Spec testReadPlink = describe "readPlink" $ do it "should read the correct Plink files" $ do @@ -172,3 +196,10 @@ let es = EigenstratIndEntry "SAMPLE0" Female "Pop1" fam = PlinkFamEntry "Pop1" "SAMPLE0" "0" "0" Female "Pop1" eigenstratInd2PlinkFam PlinkPopNameAsBoth es `shouldBe` fam + +testWriteBimInvalidId :: Spec +testWriteBimInvalidId = describe "writeBim" $ + it "should throw on SNP IDs with whitespace" $ do + let badSnps = [EigenstratSnpEntry (Chrom "11") 0 0.0 "rs0000" 'A' 'C', EigenstratSnpEntry (Chrom "11") 100000 0.001 "rs 1111" 'A' 'G'] + runSafeT (runEffect (each badSnps >-> writeBim "/tmp/invalidIdTest.bim")) `shouldThrow` + (\(SeqFormatException msg) -> "invalid SNP ID \"rs 1111\" at 11:100000" `isPrefixOf` msg) diff --git a/test/SequenceFormats/UtilsSpec.hs b/test/SequenceFormats/UtilsSpec.hs --- a/test/SequenceFormats/UtilsSpec.hs +++ b/test/SequenceFormats/UtilsSpec.hs @@ -1,13 +1,20 @@ {-# LANGUAGE OverloadedStrings #-} module SequenceFormats.UtilsSpec (spec) where -import SequenceFormats.Utils (Chrom (..), SeqFormatException (..)) +import SequenceFormats.Utils (Chrom (..), SeqFormatException (..), + decompressMultiMember) import Control.Exception (evaluate) +import qualified Data.ByteString.Char8 as B +import Pipes (each, yield) +import Pipes.GZip (compress, defaultCompression) +import qualified Pipes.Prelude as P import Test.Hspec spec :: Spec -spec = testChrom +spec = do + testChrom + testDecompressMultiMember testChrom :: Spec testChrom = describe "Chrom" $ do @@ -33,3 +40,26 @@ evaluate (Chrom "chrSSS" < Chrom "chrMT") `shouldThrow` (==SeqFormatException "cannot parse chromosome SSS") + +testDecompressMultiMember :: Spec +testDecompressMultiMember = describe "decompressMultiMember" $ do + let gz bs = P.fold (<>) B.empty id (compress defaultCompression (yield bs)) + expected = "first line\nsecond line\nthird line\n" + members <- runIO $ mapM gz ["first line\n", "second line\n", "third line\n"] + it "decompresses all members if chunks align with member boundaries" $ do + out <- P.fold (<>) B.empty id (decompressMultiMember (each members)) + out `shouldBe` expected + it "decompresses all members if chunks do not align with member boundaries" $ do + let (a, b) = B.splitAt 25 (B.concat members) + out <- P.fold (<>) B.empty id (decompressMultiMember (each [a, b])) + out `shouldBe` expected + it "throws on input that is not gzip-compressed" $ + P.fold (<>) B.empty id (decompressMultiMember (yield "this is not gzip\n")) `shouldThrow` + (== SeqFormatException "could not decompress gzip stream (zlib error code -3). The file seems to be corrupt or not gzip-compressed") + it "throws on a truncated gzip member" $ do + let truncated = B.take (B.length (B.concat members) - 5) (B.concat members) + P.fold (<>) B.empty id (decompressMultiMember (yield truncated)) `shouldThrow` + (== SeqFormatException "gzip stream ended unexpectedly. The file seems to be truncated") + it "throws on empty input" $ + P.fold (<>) B.empty id (decompressMultiMember (yield B.empty)) `shouldThrow` + (== SeqFormatException "gzip stream is empty, which is not valid gzip") diff --git a/test/SequenceFormats/VCFSpec.hs b/test/SequenceFormats/VCFSpec.hs --- a/test/SequenceFormats/VCFSpec.hs +++ b/test/SequenceFormats/VCFSpec.hs @@ -24,6 +24,7 @@ testParseVCFheader testReadVCFfromFile testReadVCFfromFileCompressed + testReadVCFfromFileMultiMember testGetGenotypes testGetDosages testIsTransversionSnp @@ -69,6 +70,19 @@ vcfRows !! 0 `shouldBe` vcf1 vcfRows !! 6 `shouldBe` vcf7 +testReadVCFfromFileMultiMember :: Spec +testReadVCFfromFileMultiMember = describe "readVCFfromFile with multi-member gzip (as in BGZF)" $ do + (vcfH, vcfRows) <- runIO . runSafeT $ do + (vcfH_, vcfProd_) <- readVCFfromFile "testDat/example.multimember.vcf.gz" + vcfRows_ <- purely P.fold list vcfProd_ + return (vcfH_, vcfRows_) + it "reads the correct sample names" $ + vcfSampleNames vcfH `shouldBe` ["12880A", "12881A", "12883A", "12884A", "12885A"] + it "reads all vcf genotype rows" $ do + length vcfRows `shouldBe` 7 + vcfRows !! 0 `shouldBe` vcf1 + vcfRows !! 6 `shouldBe` vcf7 + vcf1 :: VCFentry vcf1 = let gfields = Just (["GT", "PL"], [["0/0", "0,3,37"], ["0/0", "0,6,67"], ["0/1", "51,0,28"], ["0/0", "0,54,255"], ["0/0", "0,9,83"]]) @@ -112,10 +126,23 @@ isTransversionSnp "C" ["G"] `shouldBe` True testVcfToFreqsumEntry :: Spec -testVcfToFreqsumEntry = describe "vcfToFreqsumEntry" $ - it "should convert correctly" $ do +testVcfToFreqsumEntry = describe "vcfToFreqSumEntry" $ do + it "should convert biallelic SNPs" $ do let r = FreqSumEntry (Chrom "1") 10492 (Just "testId") Nothing 'C' 'T' [Just (0, 2), Just (0, 2), Just (1, 2), Just (0, 2), Just (0, 2)] - vcfToFreqSumEntry vcf1 `shouldReturn` r + vcfToFreqSumEntry vcf1 `shouldReturn` Just r + it "should convert sites without alternative allele and keep missing SNP IDs as Nothing" $ do + let r = FreqSumEntry (Chrom "2") 30923 Nothing Nothing 'G' 'N' [Just (2, 2), Just (2, 2), Just (2, 2), Just (2, 2), Just (2, 2)] + vcfToFreqSumEntry vcf7 `shouldReturn` Just r + it "should skip deletions" $ + vcfToFreqSumEntry vcf1 {vcfRef = "CT"} `shouldReturn` Nothing + it "should skip insertions" $ + vcfToFreqSumEntry vcf1 {vcfAlt = ["CA"]} `shouldReturn` Nothing + it "should skip multi-allelic sites" $ + vcfToFreqSumEntry vcf1 {vcfAlt = ["T", "G"]} `shouldReturn` Nothing + it "should skip spanning deletion alleles" $ + vcfToFreqSumEntry vcf1 {vcfAlt = ["*"]} `shouldReturn` Nothing + it "should still throw if genotypes are missing" $ + vcfToFreqSumEntry vcf1bad `shouldThrow` (== SeqFormatException "GT format field not found") testIsBiallelicSnp :: Spec testIsBiallelicSnp = describe "isBiallelicSnp" $ do diff --git a/testDat/example.gzipped.plink.bed b/testDat/example.gzipped.plink.bed new file mode 100644 Binary files /dev/null and b/testDat/example.gzipped.plink.bed differ diff --git a/testDat/example.multimember.vcf.gz b/testDat/example.multimember.vcf.gz new file mode 100644 Binary files /dev/null and b/testDat/example.multimember.vcf.gz differ diff --git a/testDat/example.short.plink.bed b/testDat/example.short.plink.bed new file mode 100644 --- /dev/null +++ b/testDat/example.short.plink.bed @@ -0,0 +1,1 @@ +l diff --git a/testDat/example.truncated.plink.bed b/testDat/example.truncated.plink.bed new file mode 100644 --- /dev/null +++ b/testDat/example.truncated.plink.bed @@ -0,0 +1,1 @@ +l꫋ diff --git a/testDat/example.truncated.plink.bed.gz b/testDat/example.truncated.plink.bed.gz new file mode 100644 Binary files /dev/null and b/testDat/example.truncated.plink.bed.gz differ diff --git a/testDat/example.wrongmagic.plink.bed b/testDat/example.wrongmagic.plink.bed new file mode 100644 --- /dev/null +++ b/testDat/example.wrongmagic.plink.bed @@ -0,0 +1,1 @@ +this is not a bed file