diff --git a/Changelog.md b/Changelog.md --- a/Changelog.md +++ b/Changelog.md @@ -1,5 +1,6 @@ # Changelog +- 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. - V 1.8.0.1: Allow reading arbitrary letters as reference base in Pileup format. Before only ACTG, N and M were allowed. 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.9.0.0 +version: 1.10.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 diff --git a/src/SequenceFormats/Eigenstrat.hs b/src/SequenceFormats/Eigenstrat.hs --- a/src/SequenceFormats/Eigenstrat.hs +++ b/src/SequenceFormats/Eigenstrat.hs @@ -13,27 +13,30 @@ import SequenceFormats.Utils (Chrom (..), SeqFormatException (..), consumeProducer, + deflateFinaliser, + gzipConsumer, readFileProdCheckCompress, - word, gzipConsumer, writeFromPopper) + word) import Control.Applicative ((<|>)) import Control.Exception (throw) import Control.Monad (forM_, void) import Control.Monad.Catch (MonadThrow) import Control.Monad.IO.Class (MonadIO, liftIO) +import Control.Monad.Trans.Class (lift) import qualified Data.Attoparsec.ByteString.Char8 as A import qualified Data.ByteString.Char8 as B import Data.List (isSuffixOf) import qualified Data.Streaming.Zlib as Z import Data.Vector (Vector, fromList, toList) import Pipes (Consumer, Pipe, Producer, - cat, for, yield, (>->), runEffect) + cat, for, yield, (>->)) import qualified Pipes.ByteString as PB import qualified Pipes.Prelude as P -import Pipes.Safe (MonadSafe) +import Pipes.Safe (MonadSafe, register) import qualified Pipes.Safe.Prelude as PS -import System.IO (Handle, IOMode (..), - hPutStrLn, withFile) +import System.IO (IOMode (..), hPutStrLn, + withFile) -- |A datatype to represent a single genomic SNP. The constructor arguments are: -- Chromosome, Position, Reference Allele, Alternative Allele. @@ -164,31 +167,37 @@ Unknown -> "U" -- |Function to write an Eigenstrat Snp File. Returns a consumer expecting EigenstratSnpEntries. -writeEigenstratSnp :: (MonadIO m) => Maybe Z.Deflate -- ^If Nothing, then no compression - -> Handle -- ^The Eigenstrat Snp File Handle. +writeEigenstratSnp :: (MonadSafe m) => FilePath -- ^The Eigenstrat Snp File. -> Consumer EigenstratSnpEntry m () -- ^A consumer to read EigenstratSnpEntries -writeEigenstratSnp maybeDeflate snpFileH = - let snpOutTextConsumer = case maybeDeflate of - Nothing -> PB.toHandle snpFileH - Just def -> gzipConsumer def snpFileH - toTextPipe = P.map (\(EigenstratSnpEntry chrom pos gpos gid ref alt) -> +writeEigenstratSnp snpFile = do + (_, snpFileH) <- lift $ PS.openFile snpFile WriteMode + snpOutTextConsumer <- if ".gz" `isSuffixOf` snpFile then do + def <- liftIO $ Z.initDeflate 6 (Z.WindowBits 31) + _ <- register (deflateFinaliser def snpFileH) + return $ gzipConsumer def snpFileH + else + return $ PB.toHandle snpFileH + let toTextPipe = P.map (\(EigenstratSnpEntry chrom pos gpos gid ref alt) -> 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") - in toTextPipe >-> snpOutTextConsumer + toTextPipe >-> snpOutTextConsumer -- |Function to write an Eigentrat Geno File. Returns a consumer expecting Eigenstrat Genolines. -writeEigenstratGeno :: (MonadIO m) => Maybe Z.Deflate -- ^If Nothing, then no compression - -> Handle -- ^The Genotype file handle +writeEigenstratGeno :: (MonadSafe m) => FilePath -- ^The Genotype file -> Consumer GenoLine m () -- ^A consumer to read Genotype entries. -writeEigenstratGeno maybeDeflate genoFileH = - let genoOutTextConsumer = case maybeDeflate of - Nothing -> PB.toHandle genoFileH - Just def -> gzipConsumer def genoFileH - toTextPipe = P.map (\genoLine -> +writeEigenstratGeno genoFile = do + (_, genoFileH) <- lift $ PS.openFile genoFile WriteMode + genoOutTextConsumer <- if ".gz" `isSuffixOf` genoFile then do + def <- liftIO $ Z.initDeflate 6 (Z.WindowBits 31) + _ <- register (deflateFinaliser def genoFileH) + return $ gzipConsumer def genoFileH + else + return $ PB.toHandle genoFileH + let toTextPipe = P.map (\genoLine -> let genoLineStr = B.concat . map (B.pack . show . toEigenStratNum) . toList $ genoLine in genoLineStr <> "\n") - in toTextPipe >-> genoOutTextConsumer + toTextPipe >-> genoOutTextConsumer where toEigenStratNum c = case c of HomRef -> 2 :: Int @@ -201,25 +210,10 @@ -> FilePath -- ^The Snp File -> FilePath -- ^The Ind file -> [EigenstratIndEntry] -- ^The list of individual entries - -> Producer (EigenstratSnpEntry, GenoLine) m () -- ^A consumer to read joint Snp/Genotype entries. - -> m () -writeEigenstrat genoFile snpFile indFile indEntries prod = do - snpDeflate <- if ".gz" `isSuffixOf` snpFile then fmap Just . liftIO $ Z.initDeflate 6 (Z.WindowBits 31) else return Nothing - genoDeflate <- if ".gz" `isSuffixOf` genoFile then fmap Just . liftIO $ Z.initDeflate 6 (Z.WindowBits 31) else return Nothing + -> Consumer (EigenstratSnpEntry, GenoLine) m () -- ^A consumer to read joint Snp/Genotype entries. +writeEigenstrat genoFile snpFile indFile indEntries = do liftIO $ writeEigenstratIndFile indFile indEntries - (_, snpFileH) <- PS.openFile snpFile WriteMode - (_, genoFileH) <- PS.openFile genoFile WriteMode - let snpOutConsumer = writeEigenstratSnp snpDeflate snpFileH - genoOutConsumer = writeEigenstratGeno genoDeflate genoFileH - runEffect $ prod >-> P.tee (P.map fst >-> snpOutConsumer) >-> P.map snd >-> genoOutConsumer - case snpDeflate of - Nothing -> return () - Just def -> do - let finalPop = Z.finishDeflate def - writeFromPopper finalPop snpFileH - case genoDeflate of - Nothing -> return () - Just def -> do - let finalPop = Z.finishDeflate def - writeFromPopper finalPop genoFileH - + let snpOutConsumer = writeEigenstratSnp snpFile + genoOutConsumer = writeEigenstratGeno genoFile + P.tee (P.map fst >-> snpOutConsumer) >-> P.map snd >-> genoOutConsumer + diff --git a/src/SequenceFormats/Plink.hs b/src/SequenceFormats/Plink.hs --- a/src/SequenceFormats/Plink.hs +++ b/src/SequenceFormats/Plink.hs @@ -17,13 +17,16 @@ GenoEntry (..), GenoLine, Sex (..)) import SequenceFormats.Utils (Chrom (..), consumeProducer, + deflateFinaliser, + gzipConsumer, readFileProdCheckCompress, - word, gzipConsumer, writeFromPopper) + word, writeFromPopper) import Control.Applicative ((<|>)) 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 as AB import qualified Data.Attoparsec.ByteString.Char8 as A @@ -34,14 +37,14 @@ import qualified Data.Streaming.Zlib as Z import Data.Vector (fromList, toList) import Data.Word (Word8) -import Pipes (Consumer, Producer, (>->), runEffect) +import Pipes (Consumer, Producer, (>->)) import Pipes.Attoparsec (ParsingError (..), parse) import qualified Pipes.ByteString as PB import qualified Pipes.Prelude as P -import Pipes.Safe (MonadSafe) +import Pipes.Safe (MonadSafe, register) import qualified Pipes.Safe.Prelude as PS -import System.IO (Handle, IOMode (..), - hPutStrLn, withFile) +import System.IO (IOMode (..), hPutStrLn, + withFile) -- see https://www.cog-genomics.org/plink/2.0/formats#fam data PlinkFamEntry = PlinkFamEntry { @@ -168,18 +171,21 @@ return (indEntries, P.zip snpProd genoProd) -- |Function to write a Bim file. Returns a consumer expecting EigenstratSnpEntries. -writeBim :: (MonadIO m) => Maybe Z.Deflate -- ^If Nothing, then no Compression - -> Handle -- ^The Eigenstrat Snp File handle. +writeBim :: (MonadSafe m) => FilePath -- ^The Plink Bim File. -> Consumer EigenstratSnpEntry m () -- ^A consumer to read EigenstratSnpEntries -writeBim maybeDeflate snpFileH = - let snpOutTextConsumer = case maybeDeflate of - Nothing -> PB.toHandle snpFileH - Just def -> gzipConsumer def snpFileH - toTextPipe = P.map (\(EigenstratSnpEntry chrom pos gpos gid ref alt) -> +writeBim bimFile = do + (_, bimFileH) <- lift $ PS.openFile bimFile WriteMode + bimOutTextConsumer <- if ".gz" `isSuffixOf` bimFile then do + def <- liftIO $ Z.initDeflate 6 (Z.WindowBits 31) + _ <- register (deflateFinaliser def bimFileH) + return $ gzipConsumer def bimFileH + else + return $ PB.toHandle bimFileH + let toTextPipe = P.map (\(EigenstratSnpEntry chrom pos gpos gid ref alt) -> 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") - in toTextPipe >-> snpOutTextConsumer + toTextPipe >-> bimOutTextConsumer -- |Function to write a Plink Fam file. writeFam :: (MonadIO m) => FilePath -> [PlinkFamEntry] -> m () @@ -193,20 +199,21 @@ Female -> "2" Unknown -> "0" --- |Function to write an Eigentrat Geno File. Returns a consumer expecting Eigenstrat Genolines. -writeBed :: (MonadIO m) => Maybe Z.Deflate -- ^If Nothing, then no Compression - -> Handle -- ^The Bed file handle +-- |Function to write a Plink Bed File. Returns a consumer expecting Eigenstrat Genolines. +writeBed :: (MonadSafe m) => FilePath -- ^The Bed file handle -> Consumer GenoLine m () -- ^A consumer to read Genotype entries. -writeBed maybeDeflate bedFileH = do +writeBed bedFile = do + (_, bedFileH) <- lift $ PS.openFile bedFile WriteMode let stickyBytes = BB.pack [0b01101100, 0b00011011, 0b00000001] - bedOutConsumer <- case maybeDeflate of - Nothing -> do - liftIO $ BB.hPut bedFileH stickyBytes - return $ PB.toHandle bedFileH - Just def -> do + bedOutConsumer <- if ".gz" `isSuffixOf` bedFile then do + def <- liftIO $ Z.initDeflate 6 (Z.WindowBits 31) + _ <- register (deflateFinaliser def bedFileH) pop <- liftIO (Z.feedDeflate def stickyBytes) liftIO (writeFromPopper pop bedFileH) return $ gzipConsumer def bedFileH + else do + liftIO $ BB.hPut bedFileH stickyBytes + return $ PB.toHandle bedFileH let toPlinkPipe = P.map (BB.pack . genoLineToBytes) toPlinkPipe >-> bedOutConsumer where @@ -233,24 +240,9 @@ -> FilePath -- ^The Bim File -> FilePath -- ^The Fam file -> [PlinkFamEntry] -- ^The list of individual entries - -> Producer (EigenstratSnpEntry, GenoLine) m () -- ^A consumer to read joint Snp/Genotype entries. - -> m () -writePlink bedFile bimFile famFile indEntries prod = do - bimDeflate <- if ".gz" `isSuffixOf` bimFile then fmap Just . liftIO $ Z.initDeflate 6 (Z.WindowBits 31) else return Nothing - bedDeflate <- if ".gz" `isSuffixOf` bedFile then fmap Just . liftIO $ Z.initDeflate 6 (Z.WindowBits 31) else return Nothing + -> Consumer (EigenstratSnpEntry, GenoLine) m () -- ^A consumer to read joint Snp/Genotype entries. +writePlink bedFile bimFile famFile indEntries = do liftIO $ writeFam famFile indEntries - (_, bimFileH) <- PS.openFile bimFile WriteMode - (_, bedFileH) <- PS.openFile bedFile WriteMode - let bimOutConsumer = writeBim bimDeflate bimFileH - bedOutConsumer = writeBed bedDeflate bedFileH - runEffect $ prod >-> P.tee (P.map fst >-> bimOutConsumer) >-> P.map snd >-> bedOutConsumer - case bimDeflate of - Nothing -> return () - Just def -> do - let finalPop = Z.finishDeflate def - writeFromPopper finalPop bimFileH - case bedDeflate of - Nothing -> return () - Just def -> do - let finalPop = Z.finishDeflate def - writeFromPopper finalPop bedFileH + let bimOutConsumer = writeBim bimFile + bedOutConsumer = writeBed bedFile + P.tee (P.map fst >-> bimOutConsumer) >-> P.map snd >-> bedOutConsumer diff --git a/src/SequenceFormats/Utils.hs b/src/SequenceFormats/Utils.hs --- a/src/SequenceFormats/Utils.hs +++ b/src/SequenceFormats/Utils.hs @@ -4,7 +4,7 @@ module SequenceFormats.Utils (liftParsingErrors, consumeProducer, readFileProd, readFileProdCheckCompress, - SeqFormatException(..), + SeqFormatException(..), deflateFinaliser, Chrom(..), word, gzipConsumer, writeFromPopper, Z.Deflate) where import Control.Error (readErr) @@ -17,7 +17,8 @@ import Data.Char (isSpace) import Data.List (isSuffixOf) import qualified Data.Streaming.Zlib as Z -import Pipes (Consumer, Producer, next, await) +import Pipes (Consumer, Producer, await, + next) import Pipes.Attoparsec (ParsingError (..), parsed) import qualified Pipes.ByteString as PB import Pipes.GZip (decompress) @@ -103,3 +104,8 @@ Z.PRNext bs -> do liftIO $ B.hPut h bs writeFromPopper pop h + +deflateFinaliser :: (MonadIO m) => Z.Deflate -> Handle -> m () +deflateFinaliser def h = do + let finalPop = liftIO $ Z.finishDeflate def + writeFromPopper finalPop h diff --git a/test/SequenceFormats/EigenstratSpec.hs b/test/SequenceFormats/EigenstratSpec.hs --- a/test/SequenceFormats/EigenstratSpec.hs +++ b/test/SequenceFormats/EigenstratSpec.hs @@ -4,7 +4,7 @@ import Control.Foldl (list, purely) import Control.Monad.IO.Class (liftIO) import Data.Vector (fromList) -import Pipes (each) +import Pipes (each, runEffect, (>->)) import qualified Pipes.Prelude as P import Pipes.Safe (runSafeT) import SequenceFormats.Eigenstrat (EigenstratIndEntry (..), @@ -83,7 +83,8 @@ testDatSnpProd = each mockDatEigenstratSnp testDatGenoProd = each mockDatEigenstratGeno testDatJointProd = P.zip testDatSnpProd testDatGenoProd - runSafeT $ writeEigenstrat tmpGeno tmpSnp tmpInd mockDatEigenstratInd testDatJointProd + cons = writeEigenstrat tmpGeno tmpSnp tmpInd mockDatEigenstratInd + runSafeT . runEffect $ testDatJointProd >-> cons (indEntries, esProd) <- liftIO . runSafeT $ readEigenstrat tmpGeno tmpSnp tmpInd indEntries `shouldBe` mockDatEigenstratInd snpGenoEntries <- liftIO . runSafeT $ purely P.fold list esProd @@ -99,7 +100,8 @@ testDatSnpProd = each mockDatEigenstratSnp testDatGenoProd = each mockDatEigenstratGeno testDatJointProd = P.zip testDatSnpProd testDatGenoProd - runSafeT $ writeEigenstrat tmpGeno tmpSnp tmpInd mockDatEigenstratInd testDatJointProd + cons = writeEigenstrat tmpGeno tmpSnp tmpInd mockDatEigenstratInd + runSafeT . runEffect $ testDatJointProd >-> cons (indEntries, esProd) <- liftIO . runSafeT $ readEigenstrat tmpGeno tmpSnp tmpInd indEntries `shouldBe` mockDatEigenstratInd snpGenoEntries <- liftIO . runSafeT $ purely P.fold list esProd diff --git a/test/SequenceFormats/PlinkSpec.hs b/test/SequenceFormats/PlinkSpec.hs --- a/test/SequenceFormats/PlinkSpec.hs +++ b/test/SequenceFormats/PlinkSpec.hs @@ -16,7 +16,7 @@ import Control.Foldl (list, purely) import Control.Monad.IO.Class (liftIO) import Data.Vector (fromList) -import Pipes (each) +import Pipes (each, runEffect, (>->)) import qualified Pipes.Prelude as P import Pipes.Safe (runSafeT) import Test.Hspec @@ -118,7 +118,8 @@ testDatSnpProd = each mockDatEigenstratSnp testDatGenoProd = each mockDatPlinkBed testDatJointProd = P.zip testDatSnpProd testDatGenoProd - runSafeT $ writePlink tmpGeno tmpSnp tmpInd mockDatPlinkFam testDatJointProd + cons = writePlink tmpGeno tmpSnp tmpInd mockDatPlinkFam + runSafeT . runEffect $ testDatJointProd >-> cons (indEntries, esProd) <- liftIO . runSafeT $ readPlink tmpGeno tmpSnp tmpInd indEntries `shouldBe` mockDatPlinkFam snpGenoEntries <- liftIO . runSafeT $ purely P.fold list esProd @@ -134,7 +135,8 @@ testDatSnpProd = each mockDatEigenstratSnp testDatGenoProd = each mockDatPlinkBed testDatJointProd = P.zip testDatSnpProd testDatGenoProd - runSafeT $ writePlink tmpBed tmpBim tmpFam mockDatPlinkFam testDatJointProd + cons = writePlink tmpBed tmpBim tmpFam mockDatPlinkFam + runSafeT . runEffect $ testDatJointProd >-> cons (indEntries, esProd) <- liftIO . runSafeT $ readPlink tmpBed tmpBim tmpFam indEntries `shouldBe` mockDatPlinkFam snpGenoEntries <- liftIO . runSafeT $ purely P.fold list esProd