diff --git a/Changelog.md b/Changelog.md --- a/Changelog.md +++ b/Changelog.md @@ -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. diff --git a/sequence-formats.cabal b/sequence-formats.cabal --- a/sequence-formats.cabal +++ b/sequence-formats.cabal @@ -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, diff --git a/src/SequenceFormats/Eigenstrat.hs b/src/SequenceFormats/Eigenstrat.hs --- a/src/SequenceFormats/Eigenstrat.hs +++ b/src/SequenceFormats/Eigenstrat.hs @@ -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" diff --git a/src/SequenceFormats/FreqSum.hs b/src/SequenceFormats/FreqSum.hs --- a/src/SequenceFormats/FreqSum.hs +++ b/src/SequenceFormats/FreqSum.hs @@ -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") diff --git a/src/SequenceFormats/Genomic.hs b/src/SequenceFormats/Genomic.hs --- a/src/SequenceFormats/Genomic.hs +++ b/src/SequenceFormats/Genomic.hs @@ -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 diff --git a/src/SequenceFormats/Plink.hs b/src/SequenceFormats/Plink.hs --- a/src/SequenceFormats/Plink.hs +++ b/src/SequenceFormats/Plink.hs @@ -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" diff --git a/src/SequenceFormats/VCF.hs b/src/SequenceFormats/VCF.hs --- a/src/SequenceFormats/VCF.hs +++ b/src/SequenceFormats/VCF.hs @@ -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 diff --git a/test/SequenceFormats/FreqSumSpec.hs b/test/SequenceFormats/FreqSumSpec.hs --- a/test/SequenceFormats/FreqSumSpec.hs +++ b/test/SequenceFormats/FreqSumSpec.hs @@ -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]] diff --git a/test/SequenceFormats/VCFSpec.hs b/test/SequenceFormats/VCFSpec.hs --- a/test/SequenceFormats/VCFSpec.hs +++ b/test/SequenceFormats/VCFSpec.hs @@ -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=", + "##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=", + "##contig=", + "##contig=", + "##contig=", + "##contig=", + "##contig=", + "##ALT=", + "##INFO=", + "##FORMAT=", + "##FORMAT=", + "##INFO=", + "##INFO=", + "##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 +