sequence-formats 1.10.0.0 → 1.11.0.0

raw patch · 9 files changed

+242/−112 lines, 9 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

Files

• Changelog.md +2/−0
• sequence-formats.cabal +3/−2
• src/SequenceFormats/Eigenstrat.hs +4/−4
• src/SequenceFormats/FreqSum.hs +8/−8
• src/SequenceFormats/Genomic.hs +1/−1
• src/SequenceFormats/Plink.hs +14/−14
• src/SequenceFormats/VCF.hs +113/−51
• test/SequenceFormats/FreqSumSpec.hs +7/−7
• test/SequenceFormats/VCFSpec.hs +90/−25

@@ -1,5 +1,7 @@ # Changelog +- V 1.11.0.0: Added support for writing of VCF files, including gzipping. Made some breaking API changes on top, for example+  to make the FreqSum data representation safer with respect to Ploidy. Also replaced String types in Eigenstrat and Plink formats to ByteStrings for efficiency. We anyway don't support Unicode with the AttoParsec library. - V 1.10.0.0: Brought gzip-writing support for Eigenstrat and Plink files back to non-breaking API in `writeEigenstrat` and `writePlink`. Client code can safely update from 1.8.X to 1.10.0.0. - V 1.9.0.0: Added gzip-writing support for Eigenstrat and Plink files. This required a breaking change in `writeEigenstrat` and `writePlink`. - V 1.8.1.0: Added gzip-support (read-only for now) for Plink (bed and bim files) and VCF.

@@ -1,5 +1,5 @@ name:                sequence-formats-version:             1.10.0.0+version:             1.11.0.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@@ -56,7 +56,8 @@   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, pipes-zlib+                       transformers, tasty-hunit, bytestring, containers, hspec, pipes-zlib,+                       attoparsec   other-modules:       SequenceFormats.EigenstratSpec,                        SequenceFormats.BedSpec,                        SequenceFormats.FastaSpec,

@@ -35,7 +35,7 @@ import qualified Pipes.Prelude                    as P import           Pipes.Safe                       (MonadSafe, register) import qualified Pipes.Safe.Prelude               as PS-import           System.IO                        (IOMode (..), hPutStrLn,+import           System.IO                        (IOMode (..),                                                    withFile)  -- |A datatype to represent a single genomic SNP. The constructor arguments are:@@ -52,7 +52,7 @@  -- |A datatype to represent a single individual. The constructor arguments are: -- Name, Sex and Population Name-data EigenstratIndEntry = EigenstratIndEntry String Sex String+data EigenstratIndEntry = EigenstratIndEntry B.ByteString Sex B.ByteString     deriving (Eq, Show)  -- |A datatype to represent Sex in an Eigenstrat Individual file@@ -91,7 +91,7 @@     A.skipMany1 A.space     popName <- word     void A.endOfLine-    return $ EigenstratIndEntry (B.unpack name) sex (B.unpack popName)+    return $ EigenstratIndEntry name sex popName  parseSex :: A.Parser Sex parseSex = parseMale <|> parseFemale <|> parseUnknown@@ -159,7 +159,7 @@ writeEigenstratIndFile f indEntries =     liftIO . withFile f WriteMode $ \h ->         forM_ indEntries $ \(EigenstratIndEntry name sex popName) ->-            hPutStrLn h $ name <> "\t" <> sexToStr sex <> "\t" <> popName+            B.hPutStrLn h $ name <> "\t" <> sexToStr sex <> "\t" <> popName   where     sexToStr sex = case sex of         Male    -> "M"

@@ -46,7 +46,7 @@     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.+    fsCounts     :: [Maybe (Int, Int)] -- ^A list of tuples with non-reference allele counts in each group and the total allele count. Nothing denotes missing data. } deriving (Eq, Show)  -- |This function converts a single freqSum entry to a printable freqSum line.@@ -54,9 +54,7 @@ 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-    convertToNum (Just a) = a+    countStr = B.intercalate "\t" . map (B.pack . show . maybe (-1) fst) $ maybeCounts  readFreqSumProd :: (MonadThrow m) =>     Producer B.ByteString m () -> m (FreqSumHeader, Producer FreqSumEntry m ())@@ -66,7 +64,7 @@         Nothing        -> throwM $ ParsingError [] "freqSum file exhausted"         Just (Left e)  -> throwM e         Just (Right h) -> return h-    return (header, consumeProducer parseFreqSumEntry rest)+    return (header, consumeProducer (parseFreqSumEntry (fshCounts header)) rest)  -- |A function to read a freqsum file from StdIn. Returns a pair of a freqSum Header and a Producer over all lines. readFreqSumStdIn :: (MonadIO m, MonadThrow m) => m (FreqSumHeader, Producer FreqSumEntry m ())@@ -85,11 +83,13 @@   where     tuple = (,) <$> A.takeWhile (\c -> isAlphaNum c || c == '_' || c == '-') <* A.char '(' <*> A.decimal <* A.char ')' -parseFreqSumEntry :: A.Parser FreqSumEntry-parseFreqSumEntry = FreqSumEntry <$> (Chrom <$> A.takeTill isSpace) <* A.skipSpace <*> A.decimal <*+parseFreqSumEntry :: [Int] -> A.Parser FreqSumEntry+parseFreqSumEntry denom = FreqSumEntry <$> (Chrom <$> A.takeTill isSpace) <* A.skipSpace <*> A.decimal <*     A.skipSpace <*> pure Nothing <*> pure Nothing <*> base <* A.skipSpace <*> baseOrDot <* A.skipSpace <*> counts <* A.endOfLine   where-    counts = (parseMissing <|> parseCount) `A.sepBy` A.char '\t'+    counts = do+        rawCounts <- (parseMissing <|> parseCount) `A.sepBy` A.char '\t'+        return $ map (\(c, d) -> fmap (\n -> (n, d)) c) (zip rawCounts denom)     parseMissing = A.string "-1" *> pure Nothing     parseCount = Just <$> A.decimal     base = A.satisfy (A.inClass "ACTGN")

@@ -28,7 +28,7 @@     genomicPosition (PileupRow c p _ _ _) = (c, p)  instance Genomic VCFentry where-    genomicPosition (VCFentry c p _ _ _ _ _ _ _ _) = (c, p)+    genomicPosition (VCFentry c p _ _ _ _ _ _ _) = (c, p)  chromFilter :: (Genomic e) => [Chrom] -> e -> Bool chromFilter exclusionList = (`notElem` exclusionList) . genomicChrom

@@ -33,7 +33,7 @@ import           Data.Bits                        (shiftL, shiftR, (.&.), (.|.)) import qualified Data.ByteString                  as BB import qualified Data.ByteString.Char8            as B-import           Data.List                        (intercalate, isSuffixOf)+import           Data.List                        (isSuffixOf) import qualified Data.Streaming.Zlib              as Z import           Data.Vector                      (fromList, toList) import           Data.Word                        (Word8)@@ -43,17 +43,17 @@ import qualified Pipes.Prelude                    as P import           Pipes.Safe                       (MonadSafe, register) import qualified Pipes.Safe.Prelude               as PS-import           System.IO                        (IOMode (..), hPutStrLn,+import           System.IO                        (IOMode (..),                                                    withFile)  -- see https://www.cog-genomics.org/plink/2.0/formats#fam data PlinkFamEntry = PlinkFamEntry {-    _famFamilyID     :: String,-    _famIndividualID :: String,-    _famFatherID     :: String,-    _famMotherID     :: String,+    _famFamilyID     :: B.ByteString,+    _famIndividualID :: B.ByteString,+    _famFatherID     :: B.ByteString,+    _famMotherID     :: B.ByteString,     _famSexCode      :: Sex,-    _famPhenotype    :: String+    _famPhenotype    :: B.ByteString } deriving (Eq, Show)  data PlinkPopNameMode = PlinkPopNameAsFamily | PlinkPopNameAsPhenotype | PlinkPopNameAsBoth deriving (Eq, Show)@@ -81,12 +81,12 @@ famParser :: A.Parser PlinkFamEntry famParser = do     A.skipMany A.space-    famID    <- B.unpack <$> word-    indID    <- B.unpack <$> (A.skipMany1 A.space >> word)-    fatherID <- B.unpack <$> (A.skipMany1 A.space >> word)-    motherID <- B.unpack <$> (A.skipMany1 A.space >> word)+    famID    <- word+    indID    <- A.skipMany1 A.space >> word+    fatherID <- A.skipMany1 A.space >> word+    motherID <- A.skipMany1 A.space >> word     sex      <- A.skipMany1 A.space >> parseSex-    phen     <- B.unpack <$> (A.skipMany1 A.space >> word)+    phen     <- A.skipMany1 A.space >> word     void A.endOfLine     return $ PlinkFamEntry famID indID fatherID motherID sex phen   where@@ -101,7 +101,7 @@             PlinkPopNameAsFamily    -> famId             PlinkPopNameAsPhenotype -> phen             -- If the two differ but you want both, then merge them somehow.-            PlinkPopNameAsBoth -> if famId == phen then famId else famId ++ ":" ++ phen+            PlinkPopNameAsBoth -> if famId == phen then famId else famId <> ":" <> phen     in  EigenstratIndEntry indId sex popName  eigenstratInd2PlinkFam :: PlinkPopNameMode -> EigenstratIndEntry -> PlinkFamEntry@@ -192,7 +192,7 @@ writeFam f indEntries =     liftIO . withFile f WriteMode $ \h ->         forM_ indEntries $ \(PlinkFamEntry famId indId fatherId motherId sex phen) ->-            hPutStrLn h . intercalate "\t" $ [famId, indId, fatherId, motherId, sexToStr sex, phen]+            B.hPutStrLn h . B.intercalate "\t" $ [famId, indId, fatherId, motherId, sexToStr sex, phen]   where     sexToStr sex = case sex of         Male    -> "1"

@@ -6,41 +6,52 @@  module SequenceFormats.VCF (VCFheader(..),                      VCFentry(..),+                     vcfHeaderParser,                      readVCFfromStdIn,                      readVCFfromFile,                      getGenotypes,                      getDosages,                      isTransversionSnp,                      vcfToFreqSumEntry,-                     isBiallelicSnp) where+                     isBiallelicSnp,+                     printVCFtoStdOut,+                     writeVCFfile) where  import           SequenceFormats.FreqSum          (FreqSumEntry (..)) import           SequenceFormats.Utils            (Chrom (..),                                                    SeqFormatException (..),                                                    consumeProducer,+                                                   deflateFinaliser,+                                                   gzipConsumer,                                                    readFileProdCheckCompress,-                                                   word)+                                                   word, writeFromPopper)  import           Control.Applicative              ((<|>))-import           Control.Error                    (assertErr, headErr)-import           Control.Monad                    (void)+import           Control.Error                    (atErr)+import           Control.Monad                    (forM, unless, void) import           Control.Monad.Catch              (MonadThrow, throwM)-import           Control.Monad.IO.Class           (MonadIO)+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                        (isSpace)-import           Pipes                            (Producer)+import           Data.List                        (isSuffixOf)+import           Data.Maybe                       (fromMaybe)+import qualified Data.Streaming.Zlib              as Z+import           Pipes                            (Consumer, Producer, (>->)) import           Pipes.Attoparsec                 (parse) import qualified Pipes.ByteString                 as PB-import           Pipes.Safe                       (MonadSafe)+import qualified Pipes.Prelude                    as P+import           Pipes.Safe                       (MonadSafe, register)+import qualified Pipes.Safe.Prelude               as PS+import           System.IO                        (IOMode (..))  -- |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+    vcfHeaderComments :: [B.ByteString], -- ^A list of containing all comments starting with a single '#'+    vcfSampleNames    :: [B.ByteString] -- ^The list of sample names parsed from the last comment line                              -- starting with '##'-} deriving (Show)+} deriving (Show, Eq)  -- |A Datatype representing a single VCF entry. data VCFentry = VCFentry {@@ -52,8 +63,7 @@     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.+    vcfGenotypeInfo :: Maybe ([B.ByteString], [[B.ByteString]]) -- ^An optional tuple of format tags and genotype format fields for each sample. } deriving (Show, Eq)  -- |reads a VCFheader and VCFentries from a text producer.@@ -62,7 +72,7 @@ readVCFfromProd prod = do     (res, rest) <- runStateT (parse vcfHeaderParser) prod     header <- case res of-        Nothing        -> throwM $ SeqFormatException "freqSum file exhausted"+        Nothing        -> throwM $ SeqFormatException "VCF file exhausted prematurely"         Just (Left e)  -> throwM (SeqFormatException (show e))         Just (Right h) -> return h     return (header, consumeProducer vcfEntryParser rest)@@ -76,44 +86,44 @@ readVCFfromFile = readVCFfromProd . readFileProdCheckCompress  vcfHeaderParser :: A.Parser VCFheader-vcfHeaderParser = VCFheader <$> A.many1' doubleCommentLine <*> singleCommentLine+vcfHeaderParser = VCFheader <$> A.many1' doubleCommentLine <*> (headerLineWithSamples <|> headerLineNoSamples)   where     doubleCommentLine = do         c1 <- A.string "##"         s_ <- A.takeWhile1 (/='\n')         A.endOfLine-        return . B.unpack $ c1 <> s_-    singleCommentLine = do-        void $ A.char '#'-        s_ <- A.takeWhile1 (/='\n')+        return $ c1 <> s_+    headerLineWithSamples = do+        void $ A.string "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\t"+        sampleNames <- word `A.sepBy1'` A.char '\t'         A.endOfLine-        let fields = B.splitWith (=='\t') s_-        return . drop 9 . map B.unpack $ fields+        return sampleNames+    headerLineNoSamples = A.string "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n" *> pure []  vcfEntryParser :: A.Parser VCFentry 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 <*-        A.endOfLine+        sp <*> parseInfoFields <* sp <*> parseFormatStringsAndGenotypes <* A.endOfLine     vcfEntryParserTruncated = VCFentry <$> (Chrom <$> word) <* sp <*> A.decimal <* sp <*> parseId <*         sp <*> word <* sp <*> parseAlternativeAlleles <* sp <*> parseQual <* sp <*> parseFilter <*-        sp <*> parseInfoFields <*> pure [] <*> pure [] <* A.endOfLine-    sp = A.satisfy (\c -> c == ' ' || c == '\t')+        sp <*> parseInfoFields <*> pure Nothing <* A.endOfLine+    sp = A.satisfy (\c -> c == '\t')     parseId = (parseDot *> pure Nothing) <|> (Just <$> word)     parseDot = A.char '.'     parseAlternativeAlleles = (parseDot *> pure []) <|> (parseAllele `A.sepBy1` A.char ',')-    parseAllele = A.takeTill (\c -> c == ',' || isSpace c)+    parseAllele = A.takeTill (\c -> c == ',' || c == '\t')     parseQual = (parseDot *> pure Nothing) <|> (Just <$> A.double)     parseFilter = (parseDot *> pure Nothing) <|> (Just <$> word)     parseInfoFields = (parseDot *> pure []) <|> (parseInfoField `A.sepBy1` A.char ';')-    parseInfoField = A.takeTill (\c -> c == ';' || isSpace c)+    parseInfoField = A.takeTill (\c -> c == ';' || c == '\t')+    parseFormatStringsAndGenotypes = (\f g -> Just (f, g)) <$> parseFormatStrings <* sp <*> parseGenotypeInfos      parseFormatStrings = parseFormatString `A.sepBy1` A.char ':'-    parseFormatString = A.takeTill (\c -> c == ':' || isSpace c)+    parseFormatString = A.takeTill (\c -> c == ':' || c == '\t')     parseGenotypeInfos = parseGenotype `A.sepBy1` sp     parseGenotype = parseGenoField `A.sepBy1` A.char ':'-    parseGenoField = A.takeTill (\c -> c == ':' || isSpace c)+    parseGenoField = A.takeTill (\c -> c == ':' || c == '\t' || c == '\n')  -- |returns True if the SNP is biallelic. isBiallelicSnp :: B.ByteString -> [B.ByteString] -> Bool@@ -135,34 +145,86 @@                        ((r == "C") && (a == "T")) || ((r == "T") && (a == "C"))  -- |Extracts the genotype fields (for each sapmle) from a VCF entry-getGenotypes :: VCFentry -> Either String [B.ByteString]-getGenotypes vcfEntry = do-    gtIndex <- fmap fst . headErr "GT format field not found" . filter ((=="GT") . snd) .-               zip [0..] . vcfFormatString $ vcfEntry-    return $ map (!!gtIndex) (vcfGenotypeInfo vcfEntry)+getGenotypes :: (MonadThrow m) => VCFentry -> m [B.ByteString]+getGenotypes vcfEntry = case vcfGenotypeInfo vcfEntry of+    Nothing -> throwM $ SeqFormatException "No Genotypes in this VCF"+    Just (formatField, genotypeFields) -> do+        gtIndex <- case filter ((=="GT") . snd) . zip [0..] $ formatField of+            []  -> throwM $ SeqFormatException "GT format field not found"+            [i] -> return . fst $ i+            _   -> throwM $ SeqFormatException "Multiple GT fields specified in VCF format field"+        forM genotypeFields $ \indInfo ->+            case atErr ("cannot find genotype from " ++ show indInfo) indInfo gtIndex of+                Left e  -> throwM . SeqFormatException $ e+                Right g -> return g --- |Extracts the dosages (the sum of non-reference alleles) per sample (returns a Left Error if it fails.)-getDosages :: VCFentry -> Either String [Maybe Int]+-- |Extracts the dosages (the sum of non-reference alleles) and ploidies per sample+getDosages :: (MonadThrow m) => VCFentry -> m [Maybe (Int, Int)] getDosages vcfEntry = do     genotypes <- getGenotypes vcfEntry-    let dosages = do-            gen <- genotypes-            if '.' `elem` (B.unpack gen) then-                return Nothing-            else-                return . Just $ B.count '1' gen-    return dosages+    return $ do+        gen <- genotypes+        case B.splitWith (\c -> c == '|' || c == '/') gen of+            ["0"]      -> return $ Just (0, 1)+            ["1"]      -> return $ Just (1, 1)+            ["0", "0"] -> return $ Just (0, 2)+            ["0", "1"] -> return $ Just (1, 2)+            ["1", "0"] -> return $ Just (1, 2)+            ["1", "1"] -> return $ Just (2, 2)+            _          -> return Nothing --- |Converts a VCFentry to the simpler FreqSum format (returns a Left Error if it fails.)-vcfToFreqSumEntry :: VCFentry -> Either String FreqSumEntry+-- |Converts a VCFentry to the simpler FreqSum format+vcfToFreqSumEntry :: (MonadThrow m) => VCFentry -> m FreqSumEntry vcfToFreqSumEntry vcfEntry = do-    dosages <- getDosages vcfEntry-    assertErr "multi-site reference allele" $ B.length (vcfRef vcfEntry) == 1-    assertErr "need exactly one alternative allele" $ length (vcfAlt vcfEntry) == 1-    assertErr "multi-site alternative allele" $ B.length (head . vcfAlt $ vcfEntry) == 1+    unless (B.length (vcfRef vcfEntry) == 1) . throwM $ SeqFormatException "multi-site reference allele"+    unless (length (vcfAlt vcfEntry) == 1) . throwM $ SeqFormatException "need exactly one alternative allele"+    unless (B.length (head . vcfAlt $ vcfEntry) == 1) . throwM $ SeqFormatException "multi-site alternative allele"     let ref = B.head (vcfRef vcfEntry)     let alt = B.head . head . vcfAlt $ vcfEntry-    assertErr "Invalid Reference Allele" $ ref `elem` ['A', 'C', 'T', 'G', 'N']-    assertErr "Invalid Alternative Allele" $ alt `elem` ['A', 'C', 'T', 'G', '.']+    unless (ref `elem` ['A', 'C', 'T', 'G', 'N']) . throwM $ SeqFormatException "Invalid Reference Allele"+    unless (alt `elem` ['A', 'C', 'T', 'G', '.']) . throwM $ SeqFormatException "Invalid Alternative Allele"+    dosages <- getDosages vcfEntry     return $ FreqSumEntry (vcfChrom vcfEntry) (vcfPos vcfEntry) (vcfId vcfEntry) Nothing ref alt dosages +printVCFtoStdOut :: (MonadIO m) => VCFheader -> Consumer VCFentry m ()+printVCFtoStdOut vcfh = do+    liftIO . B.putStr . vcfHeaderToText $ vcfh+    P.map vcfEntryToText >-> PB.stdout++vcfHeaderToText :: VCFheader -> B.ByteString+vcfHeaderToText (VCFheader comments names) =+    let commentsBlock = B.intercalate "\n" comments+        namesLine = case names of+            [] -> "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO"+            _  -> "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\t" <> (B.intercalate "\t" names)+    in  commentsBlock <> "\n" <> namesLine <> "\n"++vcfEntryToText :: VCFentry -> B.ByteString+vcfEntryToText e =+    let baseFieldList = [+            unChrom . vcfChrom     $ e,+            B.pack . show . vcfPos $ e,+            fromMaybe "." . vcfId  $ e,+            vcfRef e,+            if null (vcfAlt e) then "." else B.intercalate "," . vcfAlt $ e,+            maybe "." (B.pack . show) . vcfQual $ e,+            fromMaybe "." . vcfFilter $ e,+            if null (vcfInfo e) then "." else B.intercalate ";" . vcfInfo $ e]+        genotypeFieldList = case vcfGenotypeInfo e of+            Nothing -> []+            Just (f, gs) -> [B.intercalate ":" f] ++ map (B.intercalate ":") gs+    in  (<> "\n") . B.intercalate "\t" $ baseFieldList ++ genotypeFieldList++writeVCFfile :: (MonadSafe m) => FilePath -> VCFheader -> Consumer VCFentry m ()+writeVCFfile vcfFile vcfh = do+    (_, vcfFileH) <- lift $ PS.openFile vcfFile WriteMode+    vcfOutConsumer <- if ".gz" `isSuffixOf` vcfFile then do+            def <- liftIO $ Z.initDeflate 6 (Z.WindowBits 31)+            _ <- register (deflateFinaliser def vcfFileH)+            pop <- liftIO (Z.feedDeflate def (vcfHeaderToText vcfh))+            liftIO (writeFromPopper pop vcfFileH)+            return $ gzipConsumer def vcfFileH+        else do+            liftIO $ B.hPut vcfFileH (vcfHeaderToText vcfh)+            return $ PB.toHandle vcfFileH+    P.map vcfEntryToText >-> vcfOutConsumer

@@ -49,11 +49,11 @@  mockDatFsEntries :: [FreqSumEntry] mockDatFsEntries = [-    FreqSumEntry (Chrom "11") 0      Nothing Nothing 'A' 'C' [Just 1, Just 1,  Just 1, Just 1,  Just 1],-    FreqSumEntry (Chrom "11") 100000 Nothing Nothing 'A' 'G' [Just 2, Just 1,  Just 0, Just 0,  Just 0],-    FreqSumEntry (Chrom "11") 200000 Nothing Nothing 'A' 'T' [Just 0, Just 1,  Just 1, Just 1,  Just 1],-    FreqSumEntry (Chrom "11") 300000 Nothing Nothing 'C' 'A' [Just 2, Nothing, Just 1, Just 0,  Just 0],-    FreqSumEntry (Chrom "11") 400000 Nothing Nothing 'G' 'A' [Just 0, Just 1,  Just 1, Just 1,  Just 1],-    FreqSumEntry (Chrom "11") 500000 Nothing Nothing 'T' 'A' [Just 2, Just 2,  Just 1, Nothing, Just 1],-    FreqSumEntry (Chrom "11") 600000 Nothing Nothing 'G' 'T' [Just 0, Just 0,  Just 1, Nothing, Nothing]]+    FreqSumEntry (Chrom "11") 0      Nothing Nothing 'A' 'C' [Just (1, 2), Just (1, 2),  Just (1, 2), Just (1, 1),  Just (1, 1)],+    FreqSumEntry (Chrom "11") 100000 Nothing Nothing 'A' 'G' [Just (2, 2), Just (1, 2),  Just (0, 2), Just (0, 1),  Just (0, 1)],+    FreqSumEntry (Chrom "11") 200000 Nothing Nothing 'A' 'T' [Just (0, 2), Just (1, 2),  Just (1, 2), Just (1, 1),  Just (1, 1)],+    FreqSumEntry (Chrom "11") 300000 Nothing Nothing 'C' 'A' [Just (2, 2), Nothing, Just (1, 2), Just (0, 1),  Just (0, 1)],+    FreqSumEntry (Chrom "11") 400000 Nothing Nothing 'G' 'A' [Just (0, 2), Just (1, 2),  Just (1, 2), Just (1, 1),  Just (1, 1)],+    FreqSumEntry (Chrom "11") 500000 Nothing Nothing 'T' 'A' [Just (2, 2), Just (2, 2),  Just (1, 2), Nothing, Just (1, 1)],+    FreqSumEntry (Chrom "11") 600000 Nothing Nothing 'G' 'T' [Just (0, 2), Just (0, 2),  Just (1, 2), Nothing, Nothing]] 

@@ -1,19 +1,27 @@ {-# 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     (VCFentry (..), VCFheader (..),-                                          getDosages, getGenotypes,-                                          isBiallelicSnp, isTransversionSnp,-                                          readVCFfromFile, vcfToFreqSumEntry)+import           Control.Foldl                    (list, purely)+import           Data.Attoparsec.ByteString.Char8 (parseOnly)+import           Pipes                            (each, runEffect, (>->))+import qualified Pipes.Prelude                    as P+import           Pipes.Safe                       (runSafeT)+import           SequenceFormats.FreqSum          (FreqSumEntry (..))+import           SequenceFormats.Utils            (Chrom (..),+                                                   SeqFormatException (..))+import           SequenceFormats.VCF              (VCFentry (..),+                                                   VCFheader (..), getDosages,+                                                   vcfHeaderParser,+                                                   getGenotypes, isBiallelicSnp,+                                                   isTransversionSnp,+                                                   readVCFfromFile,+                                                   vcfToFreqSumEntry,+                                                   writeVCFfile) import           Test.Hspec  spec :: Spec spec = do+    testParseVCFheader     testReadVCFfromFile     testReadVCFfromFileCompressed     testGetGenotypes@@ -21,12 +29,19 @@     testIsTransversionSnp     testVcfToFreqsumEntry     testIsBiallelicSnp+    testWriteVCF +testParseVCFheader :: Spec+testParseVCFheader = describe "parseVCFheader" $ do+    let htext = "##blabla1\n##blabla2\n#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n"+    it "should correctly parse a dummy header" $ +        parseOnly vcfHeaderParser htext `shouldBe` Right (VCFheader ["##blabla1", "##blabla2"] [])+ testReadVCFfromFile :: Spec testReadVCFfromFile = describe "readVCFfromFile" $ do     (vcfH, vcfRows) <- runIO . runSafeT $ do         (vcfH_, vcfProd_) <- readVCFfromFile "testDat/example.vcf"-        vcfRows_ <- purely fold list vcfProd_+        vcfRows_ <- purely P.fold list vcfProd_         return (vcfH_, vcfRows_)     let vcfHc = vcfHeaderComments vcfH     it "reads the correct header lines" $ do@@ -42,7 +57,7 @@ testReadVCFfromFileCompressed = describe "readVCFfromFile with gzip" $ do     (vcfH, vcfRows) <- runIO . runSafeT $ do         (vcfH_, vcfProd_) <- readVCFfromFile "testDat/example.vcf.gz"-        vcfRows_ <- purely fold list vcfProd_+        vcfRows_ <- purely P.fold list vcfProd_         return (vcfH_, vcfRows_)     let vcfHc = vcfHeaderComments vcfH     it "reads the correct header lines" $ do@@ -55,29 +70,35 @@         vcfRows !! 6 `shouldBe` vcf7  vcf1 :: VCFentry-vcf1 = VCFentry (Chrom "1") 10492 (Just "testId") "C" ["T"] (Just 15.0302) Nothing ["DP=28", "PV4=1,1,0.30985,1"]-  ["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"]]+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"]])+    in  VCFentry (Chrom "1") 10492 (Just "testId") "C" ["T"] (Just 15.0302) Nothing ["DP=28", "PV4=1,1,0.30985,1"] gfields +vcf1bad :: VCFentry+vcf1bad =+    let gfields = Just (["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"]])+    in  VCFentry (Chrom "1") 10492 (Just "testId") "C" ["T"] (Just 15.0302) Nothing ["DP=28", "PV4=1,1,0.30985,1"] gfields+ vcf7 :: VCFentry-vcf7 = VCFentry (Chrom "2") 30923 Nothing "G" [] Nothing Nothing ["DP=5", "FQ=-28.9619"]-  ["GT", "PL"] [["1/1", "0,0,0"], ["1/1", "0,0,0"], ["1/1", "40,6,0"], ["1/1", "105,9,0"], ["1/1", "0,0,0"]]+vcf7 =+    let gfields = Just (["GT", "PL"], [["1/1", "0,0,0"], ["1/1", "0,0,0"], ["1/1", "40,6,0"], ["1/1", "105,9,0"], ["1/1", "0,0,0"]])+    in  VCFentry (Chrom "2") 30923 Nothing "G" [] Nothing Nothing ["DP=5", "FQ=-28.9619"] gfields  testGetGenotypes :: Spec testGetGenotypes = describe "getGenotypes" $ do     it "should successfully read genotypes if GT format field is there" $-        getGenotypes vcf1 `shouldBe` Right ["0/0", "0/0", "0/1", "0/0", "0/0"]+        getGenotypes vcf1 `shouldReturn` ["0/0", "0/0", "0/1", "0/0", "0/0"]     it "should yield Left err if GT format field isn't found" $-        getGenotypes (vcf1 {vcfFormatString=["PL"]}) `shouldBe` Left "GT format field not found"+        getGenotypes vcf1bad `shouldThrow` (== SeqFormatException "GT format field not found")  testGetDosages :: Spec testGetDosages = describe "getDosages" $ do     it "should read correct dosages" $ do-        getDosages vcf1 `shouldBe` Right [Just 0, Just 0, Just 1, Just 0, Just 0]-        let vcf1' = vcf1 {vcfGenotypeInfo=[-                ["0/0", "0,3,37"], ["0/0", "0,6,67"], [".", "51,0,28"], ["0/0", "0,54,255"],-                ["0/0", "0,9,83"]]}-        getDosages vcf1' `shouldBe` Right [Just 0, Just 0, Nothing, Just 0, Just 0]+        getDosages vcf1 `shouldReturn` [Just (0, 2), Just (0, 2), Just (1, 2), Just (0, 2), Just (0, 2)]+        let vcf1' = vcf1 {vcfGenotypeInfo = Just (["GT", "PL"], [+                ["0/0", "0,3,37"], ["0/0", "0,6,67"], [".", "51,0,28"], ["1", "0,54,255"],+                ["0/0", "0,9,83"]])}+        getDosages vcf1' `shouldReturn` [Just (0, 2), Just (0, 2), Nothing, Just (1, 1), Just (0, 2)]  testIsTransversionSnp :: Spec testIsTransversionSnp = describe "isTransversionSnp" $ do@@ -93,8 +114,8 @@ testVcfToFreqsumEntry :: Spec testVcfToFreqsumEntry = describe "vcfToFreqsumEntry" $     it "should convert correctly" $ do-        let r = Right (FreqSumEntry (Chrom "1") 10492 (Just "testId") Nothing 'C' 'T' [Just 0, Just 0, Just 1, Just 0, Just 0])-        vcfToFreqSumEntry vcf1 `shouldBe` r+        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  testIsBiallelicSnp :: Spec testIsBiallelicSnp = describe "isBiallelicSnp" $ do@@ -102,3 +123,47 @@         isBiallelicSnp "A" ["C", "T"] `shouldBe` False     it "should accept biallelic" $         isBiallelicSnp "A" ["T"] `shouldBe` True++vcfHeader :: VCFheader+vcfHeader =+    let commentLines = [+            "##fileformat=VCFv4.2",+            "##FILTER=<ID=PASS,Description=\"All filters passed\">",+            "##samtoolsVersion=1.3+htslib-1.3",+            "##samtoolsCommand=samtools mpileup -vI -f /projects1/Reference_Genomes/Human/hs37d5/hs37d5.fa -r 1:1-200000 12880A.bam 12881A.bam 12883A.bam 12884A.bam 12885A.bam",+            "##reference=file:///projects1/Reference_Genomes/Human/hs37d5/hs37d5.fa",+            "##contig=<ID=1,length=249250621>",+            "##contig=<ID=2,length=243199373>",+            "##contig=<ID=3,length=198022430>",+            "##contig=<ID=4,length=191154276>",+            "##contig=<ID=5,length=180915260>",+            "##contig=<ID=hs37d5,length=35477943>",+            "##ALT=<ID=*,Description=\"Represents allele(s) other than observed.\">",+            "##INFO=<ID=INDEL,Number=0,Type=Flag,Description=\"Indicates that the variant is an INDEL.\">",+            "##FORMAT=<ID=PL,Number=G,Type=Integer,Description=\"List of Phred-scaled genotype likelihoods\">",+            "##FORMAT=<ID=GT,Number=1,Type=String,Description=\"Genotype\">",+            "##INFO=<ID=AF1,Number=1,Type=Float,Description=\"Max-likelihood estimate of the first ALT allele frequency (assuming HWE)\">",+            "##INFO=<ID=DP4,Number=4,Type=Integer,Description=\"Number of high-quality ref-forward , ref-reverse, alt-forward and alt-reverse bases\">",+            "##bcftools_callVersion=1.3+htslib-1.3",+            "##bcftools_callCommand=call -c -v"]+        sampleNames = ["12880A", "12881A", "12883A", "12884A", "12885A"]+    in  VCFheader commentLines sampleNames++testWriteVCF :: Spec+testWriteVCF = describe "writeVCF" $ do+    let tmpVCF = "/tmp/vcfWriteTest.vcf"+        testDatVCFprod = each [vcf1, vcf7]+        cons = writeVCFfile tmpVCF vcfHeader+    runIO . runSafeT . runEffect $ testDatVCFprod >-> cons+    (vcfH, vcfRows) <- runIO . runSafeT $ do+        (vcfH_, vcfProd_) <- readVCFfromFile tmpVCF+        vcfRows_ <- purely P.fold list vcfProd_+        return (vcfH_, vcfRows_)+    it "correctly write and reads back VCF data" $ do+        let vcfHc = vcfHeaderComments vcfH+        vcfHc !! 0 `shouldBe` "##fileformat=VCFv4.2"+        vcfHc !! 18 `shouldBe` "##bcftools_callCommand=call -c -v"+        vcfSampleNames vcfH `shouldBe` ["12880A", "12881A", "12883A", "12884A", "12885A"]+        vcfRows !! 0 `shouldBe` vcf1+        vcfRows !! 1 `shouldBe` vcf7+