diff --git a/Changelog.md b/Changelog.md --- a/Changelog.md +++ b/Changelog.md @@ -1,3 +1,5 @@ +V 1.6.0: Plink Output now supported. + V 1.5.1.4: Plink BIM parsing now automatically converts from Numbers (0,1,2,3,4) to Letters (N,A,C,T,G). V 1.5.1.3: added possibility to parse allele names 01234 in bim files. 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.5.2 +version: 1.6.0 license: GPL-3 license-file: LICENSE maintainer: stephan.schiffels@mac.com @@ -47,19 +47,19 @@ default-language: Haskell2010 build-depends: base >=4.7 && <5, - containers >=0.6.0.1, + containers >=0.6.2.1, errors >=2.3.0, - attoparsec >=0.13.2.2, - pipes >=4.3.11, + attoparsec >=0.13.2.4, + pipes >=4.3.14, transformers >=0.5.6.2, - bytestring >=0.10.8.2, - lens-family >=1.2.3, + bytestring >=0.10.12.0, + lens-family >=2.0.0, pipes-bytestring >=2.1.6, - foldl >=1.4.5, - exceptions >=0.10.2, - pipes-safe >=2.3.1, + foldl >=1.4.10, + exceptions >=0.10.4, + pipes-safe >=2.3.2, pipes-attoparsec >=0.5.1.5, - vector >=0.12.0.3 + vector >=0.12.1.2 test-suite sequenceFormatTests type: exitcode-stdio-1.0 @@ -77,15 +77,15 @@ default-language: Haskell2010 build-depends: - base >=4.12.0.0, + base >=4.14.1.0, sequence-formats -any, - foldl >=1.4.5, - pipes >=4.3.11, - pipes-safe >=2.3.1, + foldl >=1.4.10, + pipes >=4.3.14, + pipes-safe >=2.3.2, tasty >=1.2.3, - vector >=0.12.0.3, + vector >=0.12.1.2, transformers >=0.5.6.2, - tasty-hunit >=0.10.0.2, - bytestring >=0.10.8.2, - containers >=0.6.0.1, - hspec >=2.7.1 + tasty-hunit >=0.10.0.3, + bytestring >=0.10.12.0, + containers >=0.6.2.1, + hspec >=2.7.8 diff --git a/src/SequenceFormats/Plink.hs b/src/SequenceFormats/Plink.hs --- a/src/SequenceFormats/Plink.hs +++ b/src/SequenceFormats/Plink.hs @@ -1,6 +1,6 @@ {-# LANGUAGE BinaryLiterals #-} {-# LANGUAGE OverloadedStrings #-} -module SequenceFormats.Plink (readBimStdIn, readBimFile, writeBim, readFamFile, readPlinkBedFile, readPlink) where +module SequenceFormats.Plink (readBimStdIn, readBimFile, writeBim, readFamFile, readPlinkBedFile, readPlink, writePlink) where import SequenceFormats.Eigenstrat (EigenstratIndEntry (..), EigenstratSnpEntry (..), @@ -10,24 +10,26 @@ readFileProd, word) import Control.Applicative ((<|>)) -import Control.Monad (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 as AB import qualified Data.Attoparsec.ByteString.Char8 as A -import Data.Bits (shiftR, (.&.)) +import Data.Bits (shiftL, shiftR, (.&.), (.|.)) import qualified Data.ByteString as BB import qualified Data.ByteString.Char8 as B -import Data.Vector (fromList) +import Data.Vector (fromList, toList) +import Data.Word (Word8) 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 qualified Pipes.Safe.Prelude +import qualified Pipes.Safe.Prelude as PS import System.IO (Handle, IOMode (..), - withFile) + hPutStrLn, withFile) bimParser :: A.Parser EigenstratSnpEntry @@ -98,7 +100,7 @@ -- |A function to read a bed file from a file. Returns a Producer over all lines. readPlinkBedFile :: (MonadSafe m) => FilePath -> Int -> m (Producer GenoLine m ()) -readPlinkBedFile file nrInds = readPlinkBedProd nrInds (Pipes.Safe.Prelude.withFile file ReadMode PB.fromHandle) +readPlinkBedFile file nrInds = readPlinkBedProd nrInds (PS.withFile file ReadMode PB.fromHandle) -- |Function to read a Bim File from StdIn. Returns a Pipes-Producer over the EigenstratSnpEntries. readBimStdIn :: (MonadThrow m, MonadIO m) => Producer EigenstratSnpEntry m () @@ -132,5 +134,58 @@ writeBim snpFileH = let snpOutTextConsumer = PB.toHandle snpFileH toTextPipe = P.map (\(EigenstratSnpEntry chrom pos gpos gid ref alt) -> - B.intercalate "\t" [unChrom chrom, gid, B.pack (show gpos), B.pack (show pos), B.singleton ref, B.singleton alt]) + 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 + +-- |Function to write a Plink Fam file. +writeFam :: (MonadIO m) => FilePath -> [EigenstratIndEntry] -> m () +writeFam f indEntries = + liftIO . withFile f WriteMode $ \h -> + forM_ indEntries $ \(EigenstratIndEntry name sex popName) -> + hPutStrLn h $ popName <> "\t" <> name <> "\t0\t0\t" <> sexToStr sex <> "\t0" + where + sexToStr sex = case sex of + Male -> "1" + Female -> "2" + Unknown -> "0" + +-- |Function to write an Eigentrat Geno File. Returns a consumer expecting Eigenstrat Genolines. +writeBed :: (MonadIO m) => Handle -- ^The Bed file handle + -> Consumer GenoLine m () -- ^A consumer to read Genotype entries. +writeBed bedFileH = do + liftIO $ BB.hPut bedFileH (BB.pack [0b01101100, 0b00011011, 0b00000001]) + let bedOutConsumer = PB.toHandle bedFileH + toPlinkPipe = P.map (BB.pack . genoLineToBytes) + toPlinkPipe >-> bedOutConsumer + where + genoLineToBytes :: GenoLine -> [Word8] + genoLineToBytes genoLine = go (toList genoLine) + where + go :: [GenoEntry] -> [Word8] + go (g1 : g2 : g3 : g4 : rest) = constructByte [g1, g2, g3, g4] : go rest + -- go (g1 : g2 : g3 : g4 : rest) = constructByte [g4, g3, g2, g1] : go rest + go genoEntries = [constructByte genoEntries] + constructByte :: [GenoEntry] -> Word8 + constructByte [] = error "constructByte - should never happen" + constructByte [g] = genoEntryToByte g + constructByte (g:gs) = shiftL (constructByte gs) 2 .|. genoEntryToByte g + +genoEntryToByte :: GenoEntry -> Word8 +genoEntryToByte HomRef = 0b00000000 +genoEntryToByte HomAlt = 0b00000011 +genoEntryToByte Het = 0b00000010 +genoEntryToByte Missing = 0b00000001 + +-- |Function to write a Plink Database. Returns a consumer expecting joint Snp- and Genotype lines. +writePlink :: (MonadSafe m) => FilePath -- ^The Bed file + -> FilePath -- ^The Bim File + -> FilePath -- ^The Fam file + -> [EigenstratIndEntry] -- ^The list of individual entries + -> Consumer (EigenstratSnpEntry, GenoLine) m () -- ^A consumer to read joint Snp/Genotype entries. +writePlink bedFile bimFile famFile indEntries = do + liftIO $ writeFam famFile indEntries + let bimOutConsumer = PS.withFile bimFile WriteMode writeBim + bedOutConsumer = PS.withFile bedFile WriteMode writeBed + P.tee (P.map fst >-> bimOutConsumer) >-> P.map snd >-> bedOutConsumer diff --git a/test/SequenceFormats/EigenstratSpec.hs b/test/SequenceFormats/EigenstratSpec.hs --- a/test/SequenceFormats/EigenstratSpec.hs +++ b/test/SequenceFormats/EigenstratSpec.hs @@ -72,3 +72,4 @@ indEntries `shouldBe` mockDatEigenstratInd snpGenoEntries <- liftIO . runSafeT $ purely P.fold list esProd (map fst snpGenoEntries) `shouldBe` mockDatEigenstratSnp + (map snd snpGenoEntries) `shouldBe` mockDatEigenstratGeno diff --git a/test/SequenceFormats/PlinkSpec.hs b/test/SequenceFormats/PlinkSpec.hs --- a/test/SequenceFormats/PlinkSpec.hs +++ b/test/SequenceFormats/PlinkSpec.hs @@ -2,12 +2,17 @@ module SequenceFormats.PlinkSpec (spec) where import SequenceFormats.Eigenstrat (EigenstratIndEntry (..), - EigenstratSnpEntry (..), Sex (..), GenoLine, GenoEntry(..)) -import SequenceFormats.Plink (readBimFile, readFamFile, readPlinkBedFile) + EigenstratSnpEntry (..), + GenoEntry (..), GenoLine, Sex (..)) +import SequenceFormats.Plink (readBimFile, readFamFile, + readPlink, readPlinkBedFile, + writePlink) import SequenceFormats.Utils (Chrom (..)) import Control.Foldl (list, purely) +import Control.Monad.IO.Class (liftIO) import Data.Vector (fromList) +import Pipes (each, runEffect, (>->)) import qualified Pipes.Prelude as P import Pipes.Safe (runSafeT) import Test.Hspec @@ -17,6 +22,8 @@ testReadBimFile testReadFamFile testReadBedFile + testReadPlink + testWritePlink mockDatEigenstratSnp :: [EigenstratSnpEntry] mockDatEigenstratSnp = [ @@ -65,3 +72,32 @@ bedProd <- readPlinkBedFile fn 5 purely P.fold list bedProd bedDat `shouldBe` mockDatPlinkBed + +testReadPlink :: Spec +testReadPlink = describe "readPlink" $ do + it "should read the correct Plink files" $ do + let bimFile = "testDat/example.plink.bim" + famFile = "testDat/example.plink.fam" + bedFile = "testDat/example.plink.bed" + (indEntries, esProd) <- runSafeT $ readPlink bedFile bimFile famFile + indEntries `shouldBe` mockDatEigenstratInd + snpGenoEntries <- runSafeT $ purely P.fold list esProd + (map fst snpGenoEntries) `shouldBe` mockDatEigenstratSnp + (map snd snpGenoEntries) `shouldBe` mockDatPlinkBed + +testWritePlink :: Spec +testWritePlink = describe "writePlink" $ do + it "should write and read back Plink data correctly" $ do + let tmpGeno = "/tmp/plinkWriteTest.bed" + tmpSnp = "/tmp/plinkWriteTest.bim" + tmpInd = "/tmp/plinkWriteTest.fam" + testDatSnpProd = each mockDatEigenstratSnp + testDatGenoProd = each mockDatPlinkBed + testDatJointProd = P.zip testDatSnpProd testDatGenoProd + liftIO . runSafeT . runEffect $ + testDatJointProd >-> writePlink tmpGeno tmpSnp tmpInd mockDatEigenstratInd + (indEntries, esProd) <- liftIO . runSafeT $ readPlink tmpGeno tmpSnp tmpInd + indEntries `shouldBe` mockDatEigenstratInd + snpGenoEntries <- liftIO . runSafeT $ purely P.fold list esProd + (map fst snpGenoEntries) `shouldBe` mockDatEigenstratSnp + (map snd snpGenoEntries) `shouldBe` mockDatPlinkBed diff --git a/testDat/example.fam b/testDat/example.fam --- a/testDat/example.fam +++ b/testDat/example.fam @@ -1,5 +1,5 @@ - 1 SAMPLE0 0 0 2 2 - 2 SAMPLE1 0 0 1 2 - 3 SAMPLE2 0 0 2 1 - 4 SAMPLE3 0 0 1 1 - 5 SAMPLE4 0 0 2 1 + 1 SAMPLE0 0 0 2 0 + 2 SAMPLE1 0 0 1 0 + 3 SAMPLE2 0 0 2 0 + 4 SAMPLE3 0 0 1 0 + 5 SAMPLE4 0 0 2 0