sequence-formats 1.5.2 → 1.6.0

raw patch · 6 files changed

+128/−34 lines, 6 filesdep ~attoparsecdep ~basedep ~bytestringPVP ok

version bump matches the API change (PVP)

Dependency ranges changed: attoparsec, base, bytestring, containers, exceptions, foldl, hspec, lens-family, pipes, pipes-safe, tasty-hunit, vector

API changes (from Hackage documentation)

Files

• Changelog.md +2/−0
• sequence-formats.cabal +19/−19
• src/SequenceFormats/Plink.hs +63/−8
• test/SequenceFormats/EigenstratSpec.hs +1/−0
• test/SequenceFormats/PlinkSpec.hs +38/−2
• testDat/example.fam +5/−5

@@ -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.

@@ -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

@@ -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

@@ -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

@@ -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

@@ -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