sequence-formats 1.6.6.1 → 1.6.7.0

raw patch · 7 files changed

+123/−43 lines, 7 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• Changelog.md +8/−0
• sequence-formats.cabal +1/−1
• src/SequenceFormats/Pileup.hs +3/−3
• src/SequenceFormats/Plink.hs +57/−19
• test/SequenceFormats/PlinkSpec.hs +47/−13
• testDat/example.fam +5/−5
• testDat/example.pileup +2/−2

@@ -1,7 +1,15 @@+V 1.6.7.0:+* Added an option to Plink output, which allows users to store the population mame in either the first or last column, or both+* Fixed a bug in the pileup-parser which failed at reference-skip symbols (">" and "<").+ V 1.6.6.0: Added a minimal Bed parser and a bed-filter+ V 1.6.5.0: Long chromosome names are now parsed correctly+ V 1.6.4: Fixed a bug in Fasta parsing which would skip the first line if the header was minimal.+ V 1.6.3: Pileup reference-base is read as upper case, even if it's not in the file.+ V 1.6.2: Made compatible with GHC 9 (thanks to github user sjakobi)  V 1.6.0: Plink Output now supported.

@@ -1,6 +1,6 @@ cabal-version:      >=1.10 name:               sequence-formats-version:            1.6.6.1+version:            1.6.7.0 license:            GPL-3 license-file:       LICENSE maintainer:         stephan.schiffels@mac.com

@@ -68,9 +68,9 @@         | x == ',' = (refA, ReverseStrand) : go xs         | x `elem` ("ACTGN" :: String) = (x, ForwardStrand) : go xs         | x `elem` ("actgn" :: String) = (toUpper x, ReverseStrand) : go xs-        | x `elem` ("$*" :: String) = go xs-        | x == '^' = go (drop 1 xs)-        | (x == '+' || x == '-') =+        | x `elem` ("$*#<>" :: String) = go xs+        | x == '^' = go (drop 1 xs) -- skip the next character, which is the mapping quality +        | x == '+' || x == '-' =  -- insertions or deletions, followed by a decimal number             let [(num, rest)] = reads xs in go (drop num rest)         | otherwise = error $ "cannot parse read base string: " ++ (x:xs)     go [] = []

@@ -1,6 +1,16 @@ {-# LANGUAGE BinaryLiterals    #-} {-# LANGUAGE OverloadedStrings #-}-module SequenceFormats.Plink (readBimStdIn, readBimFile, writeBim, readFamFile, readPlinkBedFile, readPlink, writePlink) where+module SequenceFormats.Plink (readBimStdIn,+                              readBimFile,+                              writeBim,+                              readFamFile,+                              readPlinkBedFile,+                              readPlink,+                              writePlink,+                              PlinkFamEntry(..),+                              plinkFam2EigenstratInd,+                              eigenstratInd2PlinkFam,+                              PlinkPopNameMode(..)) where  import           SequenceFormats.Eigenstrat       (EigenstratIndEntry (..),                                                    EigenstratSnpEntry (..),@@ -19,6 +29,7 @@ import           Data.Bits                        (shiftL, shiftR, (.&.), (.|.)) import qualified Data.ByteString                  as BB import qualified Data.ByteString.Char8            as B+import           Data.List                        (intercalate) import           Data.Vector                      (fromList, toList) import           Data.Word                        (Word8) import           Pipes                            (Consumer, Producer, (>->))@@ -30,7 +41,18 @@ import           System.IO                        (Handle, IOMode (..),                                                    hPutStrLn, withFile) +-- see https://www.cog-genomics.org/plink/2.0/formats#fam+data PlinkFamEntry = PlinkFamEntry {+    _famFamilyID     :: String,+    _famIndividualID :: String,+    _famFatherID     :: String,+    _famMotherID     :: String,+    _famSexCode      :: Sex,+    _famPhenotype    :: String+} deriving (Eq, Show) +data PlinkPopNameMode = PlinkPopNameAsFamily | PlinkPopNameAsPhenotype | PlinkPopNameAsBoth deriving (Eq, Show)+ bimParser :: A.Parser EigenstratSnpEntry bimParser = do     chrom      <- word@@ -51,23 +73,39 @@     convertNum '4' = 'T'     convertNum x   = x -famParser :: A.Parser EigenstratIndEntry+famParser :: A.Parser PlinkFamEntry famParser = do     A.skipMany A.space-    pop <- word-    ind <- A.skipMany1 A.space >> word-    _   <- A.skipMany1 A.space >> A.decimal :: A.Parser Int-    _   <- A.skipMany1 A.space >> A.decimal :: A.Parser Int-    sex <- A.skipMany1 A.space >> parseSex-    _   <- A.skipMany1 A.space >> word+    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)+    sex      <- A.skipMany1 A.space >> parseSex+    phen     <- B.unpack <$> (A.skipMany1 A.space >> word)     void A.endOfLine-    return $ EigenstratIndEntry (B.unpack ind) sex (B.unpack pop)+    return $ PlinkFamEntry famID indID fatherID motherID sex phen   where     parseSex = parseMale <|> parseFemale <|> parseUnknown     parseMale = A.char '1' >> return Male     parseFemale = A.char '2' >> return Female     parseUnknown = A.anyChar >> return Unknown +plinkFam2EigenstratInd :: PlinkPopNameMode -> PlinkFamEntry -> EigenstratIndEntry+plinkFam2EigenstratInd plinkPopNameMode (PlinkFamEntry famId indId _ _ sex phen) =+    let popName = case plinkPopNameMode of+            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+    in  EigenstratIndEntry indId sex popName++eigenstratInd2PlinkFam :: PlinkPopNameMode -> EigenstratIndEntry -> PlinkFamEntry+eigenstratInd2PlinkFam plinkPopNameMode (EigenstratIndEntry indId sex popName)=+    case plinkPopNameMode of+        PlinkPopNameAsFamily    -> PlinkFamEntry popName indId "0" "0" sex "0"+        PlinkPopNameAsPhenotype -> PlinkFamEntry "DummyFamily" indId "0" "0" sex popName+        PlinkPopNameAsBoth      -> PlinkFamEntry popName indId "0" "0" sex popName+ bedHeaderParser :: AB.Parser () bedHeaderParser = do     void $ AB.word8 0b01101100 -- magic number I for BED files@@ -110,7 +148,7 @@ readBimFile = consumeProducer bimParser . readFileProd  -- |Function to read a Plink fam file. Returns the Eigenstrat Individual Entries as list.-readFamFile :: (MonadIO m) => FilePath -> m [EigenstratIndEntry]+readFamFile :: (MonadIO m) => FilePath -> m [PlinkFamEntry] readFamFile fn =     liftIO . withFile fn ReadMode $ \handle ->         P.toListM $ consumeProducer famParser (PB.fromHandle handle)@@ -119,7 +157,7 @@ readPlink :: (MonadSafe m) => FilePath -- ^The Bed file                -> FilePath -- ^The Bim File                -> FilePath -- ^The Fam file-               -> m ([EigenstratIndEntry], Producer (EigenstratSnpEntry, GenoLine) m ()) -- The return pair of individual entries and a joint Snp/Geno Producer.+               -> m ([PlinkFamEntry], Producer (EigenstratSnpEntry, GenoLine) m ()) -- The return pair of individual entries and a joint Snp/Geno Producer. readPlink bedFile bimFile famFile = do     indEntries <- readFamFile famFile     let nrInds = length indEntries@@ -139,11 +177,11 @@     in  toTextPipe >-> snpOutTextConsumer  -- |Function to write a Plink Fam file.-writeFam :: (MonadIO m) => FilePath -> [EigenstratIndEntry] -> m ()+writeFam :: (MonadIO m) => FilePath -> [PlinkFamEntry] -> 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"+        forM_ indEntries $ \(PlinkFamEntry famId indId fatherId motherId sex phen) ->+            hPutStrLn h . intercalate "\t" $ [famId, indId, fatherId, motherId, sexToStr sex, phen]   where     sexToStr sex = case sex of         Male    -> "1"@@ -163,12 +201,12 @@     genoLineToBytes genoLine = go (toList genoLine)       where         go :: [GenoEntry] -> [Word8]-        go [] = [] -- empty list for recursion stop+        go []                         = [] -- empty list for recursion stop         go (g1 : g2 : g3 : g4 : rest) = constructByte [g1, g2, g3, g4] : go rest -- at least 5 entries -> more than 1 byte-        go genoEntries = [constructByte genoEntries] -- four or less entries -> 1 byte+        go genoEntries                = [constructByte genoEntries] -- four or less entries -> 1 byte         constructByte :: [GenoEntry] -> Word8-        constructByte [] = error "constructByte - should never happen"-        constructByte [g] = genoEntryToByte g+        constructByte []     = error "constructByte - should never happen"+        constructByte [g]    = genoEntryToByte g         constructByte (g:gs) = shiftL (constructByte gs) 2 .|. genoEntryToByte g  genoEntryToByte :: GenoEntry -> Word8@@ -181,7 +219,7 @@ writePlink :: (MonadSafe m) => FilePath -- ^The Bed file                 -> FilePath -- ^The Bim File                 -> FilePath -- ^The Fam file-                -> [EigenstratIndEntry] -- ^The list of individual entries+                -> [PlinkFamEntry] -- ^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

@@ -6,7 +6,9 @@                                              GenoEntry (..), GenoLine, Sex (..)) import           SequenceFormats.Plink      (readBimFile, readFamFile,                                              readPlink, readPlinkBedFile,-                                             writePlink)+                                             writePlink, PlinkFamEntry(..),+                                             plinkFam2EigenstratInd, eigenstratInd2PlinkFam,+                                             PlinkPopNameMode(..)) import           SequenceFormats.Utils      (Chrom (..))  import           Control.Foldl              (list, purely)@@ -24,6 +26,8 @@     testReadBedFile     testReadPlink     testWritePlink+    testFam2Ind+    testInd2Fam  mockDatEigenstratSnp :: [EigenstratSnpEntry] mockDatEigenstratSnp = [@@ -35,13 +39,13 @@     EigenstratSnpEntry (Chrom "11") 500000 0.005000 "rs5555" 'T' 'A',     EigenstratSnpEntry (Chrom "11") 600000 0.006000 "rs6666" 'G' 'T'] -mockDatEigenstratInd :: [EigenstratIndEntry]-mockDatEigenstratInd = [-    EigenstratIndEntry "SAMPLE0" Female "1",-    EigenstratIndEntry "SAMPLE1" Male "2",-    EigenstratIndEntry "SAMPLE2" Female "3",-    EigenstratIndEntry "SAMPLE3" Male "4",-    EigenstratIndEntry "SAMPLE4" Female "5"]+mockDatPlinkFam :: [PlinkFamEntry]+mockDatPlinkFam = [+    PlinkFamEntry "Family1" "SAMPLE0" "0" "0" Female "Pop1",+    PlinkFamEntry "Family2" "SAMPLE1" "0" "0" Male "Pop2",+    PlinkFamEntry "Family3" "SAMPLE2" "0" "0" Female "Pop3",+    PlinkFamEntry "Family4" "SAMPLE3" "0" "0" Male "Pop4",+    PlinkFamEntry "Family5" "SAMPLE4" "0" "0" Female "Pop5"]  mockDatPlinkBed :: [GenoLine] mockDatPlinkBed = [@@ -62,7 +66,7 @@ testReadFamFile :: Spec testReadFamFile = describe "readFamFile" $     it "should read a FAM file correctly" $ do-        readFamFile "testDat/example.fam" `shouldReturn` mockDatEigenstratInd+        readFamFile "testDat/example.fam" `shouldReturn` mockDatPlinkFam  testReadBedFile :: Spec testReadBedFile = describe "readBedFile" $@@ -77,10 +81,10 @@ testReadPlink = describe "readPlink" $ do     it "should read the correct Plink files" $ do         let bimFile = "testDat/example.plink.bim"-            famFile = "testDat/example.plink.fam"+            famFile = "testDat/example.fam"             bedFile = "testDat/example.plink.bed"         (indEntries, esProd) <- runSafeT $ readPlink bedFile bimFile famFile-        indEntries `shouldBe` mockDatEigenstratInd+        indEntries `shouldBe` mockDatPlinkFam         snpGenoEntries <- runSafeT $ purely P.fold list esProd         (map fst snpGenoEntries) `shouldBe` mockDatEigenstratSnp         (map snd snpGenoEntries) `shouldBe` mockDatPlinkBed@@ -95,9 +99,39 @@             testDatGenoProd = each mockDatPlinkBed             testDatJointProd = P.zip testDatSnpProd testDatGenoProd         liftIO . runSafeT . runEffect $-            testDatJointProd >-> writePlink tmpGeno tmpSnp tmpInd mockDatEigenstratInd+            testDatJointProd >-> writePlink tmpGeno tmpSnp tmpInd mockDatPlinkFam         (indEntries, esProd) <- liftIO . runSafeT $ readPlink tmpGeno tmpSnp tmpInd-        indEntries `shouldBe` mockDatEigenstratInd+        indEntries `shouldBe` mockDatPlinkFam         snpGenoEntries <- liftIO . runSafeT $ purely P.fold list esProd         (map fst snpGenoEntries) `shouldBe` mockDatEigenstratSnp         (map snd snpGenoEntries) `shouldBe` mockDatPlinkBed++testFam2Ind :: Spec+testFam2Ind = describe "plinkFam2EigenstratInd" $ do+    it "should correctly convert with family-option" $ do+        let fam = PlinkFamEntry "Family1" "SAMPLE0" "0" "0" Female "Pop1"+        plinkFam2EigenstratInd PlinkPopNameAsFamily fam `shouldBe` EigenstratIndEntry "SAMPLE0" Female "Family1"+    it "should correctly convert with phenotype-option" $ do+        let fam = PlinkFamEntry "Family1" "SAMPLE0" "0" "0" Female "Pop1"+        plinkFam2EigenstratInd PlinkPopNameAsPhenotype fam `shouldBe` EigenstratIndEntry "SAMPLE0" Female "Pop1"+    it "should correctly convert with both-option when they differ" $ do+        let fam = PlinkFamEntry "Family1" "SAMPLE0" "0" "0" Female "Pop1"+        plinkFam2EigenstratInd PlinkPopNameAsBoth fam `shouldBe` EigenstratIndEntry "SAMPLE0" Female "Family1:Pop1"+    it "should correctly convert with both-option when they are the same" $ do+        let fam = PlinkFamEntry "Pop1" "SAMPLE0" "0" "0" Female "Pop1"+        plinkFam2EigenstratInd PlinkPopNameAsBoth fam `shouldBe` EigenstratIndEntry "SAMPLE0" Female "Pop1"+    +testInd2Fam :: Spec+testInd2Fam = describe "eigenstratInd2PlinkFam" $ do+    it "should correctly convert with family-option" $ do+        let es = EigenstratIndEntry "SAMPLE0" Female "Pop1"+            fam = PlinkFamEntry "Pop1" "SAMPLE0" "0" "0" Female "0"+        eigenstratInd2PlinkFam PlinkPopNameAsFamily es `shouldBe` fam+    it "should correctly convert with phenotype-option" $ do+        let es = EigenstratIndEntry "SAMPLE0" Female "Pop1"+            fam = PlinkFamEntry "DummyFamily" "SAMPLE0" "0" "0" Female "Pop1"+        eigenstratInd2PlinkFam PlinkPopNameAsPhenotype es `shouldBe` fam+    it "should correctly convert with both-option" $ do+        let es = EigenstratIndEntry "SAMPLE0" Female "Pop1"+            fam = PlinkFamEntry "Pop1" "SAMPLE0" "0" "0" Female "Pop1"+        eigenstratInd2PlinkFam PlinkPopNameAsBoth es `shouldBe` fam

@@ -1,5 +1,5 @@-     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+     Family1      SAMPLE0 0 0 2 Pop1+     Family2      SAMPLE1 0 0 1 Pop2+     Family3      SAMPLE2 0 0 2 Pop3+     Family4      SAMPLE3 0 0 1 Pop4+     Family5      SAMPLE4 0 0 2 Pop5

@@ -1,3 +1,3 @@ 1 1000 A 4 ..,C*. IIII 3 ..ac* III 8 ,,.,ccAA-3acccA IIIIIIII-1 2000 c 4 ..C+4IIIIa IIII 3 ACt..* III 8 .a.a.C.C* IIIIIIII-2 1000 G 4 ,,.,,..-3AccC*. IIII 3 .,,* III 8 .,,.,,* IIIIIIII+1 2000 c 4 ..C+4IIIIa< IIII 3 ACt..* III 8 .a.a.C.C* IIIIIIII+2 1000 G 4 ,,.,,.#.-3Acc>C*. IIII 3 .,,* III 8 .,,.,,* IIIIIIII