diff --git a/Changelog.md b/Changelog.md --- a/Changelog.md +++ b/Changelog.md @@ -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. 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.6.6.1 +version: 1.6.7.0 license: GPL-3 license-file: LICENSE maintainer: stephan.schiffels@mac.com diff --git a/src/SequenceFormats/Pileup.hs b/src/SequenceFormats/Pileup.hs --- a/src/SequenceFormats/Pileup.hs +++ b/src/SequenceFormats/Pileup.hs @@ -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 [] = [] 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,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 diff --git a/test/SequenceFormats/PlinkSpec.hs b/test/SequenceFormats/PlinkSpec.hs --- a/test/SequenceFormats/PlinkSpec.hs +++ b/test/SequenceFormats/PlinkSpec.hs @@ -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 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 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 diff --git a/testDat/example.pileup b/testDat/example.pileup --- a/testDat/example.pileup +++ b/testDat/example.pileup @@ -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