diff --git a/Changelog.md b/Changelog.md --- a/Changelog.md +++ b/Changelog.md @@ -26,4 +26,9 @@ V 1.3.2.1: Fixed a hard-coded absolute path in the test-suite. -V. 1.3.3: Added Pileup as new format. Changed all tests to Hspec. +V 1.3.3: Added Pileup as new format. Changed all tests to Hspec. + +V 1.4.0: Added three features: + - Chromosomes now include X, Y and MT (or chrX, chrY, chrMT), in that order after chr22. + - SNP rsId information is now internally included as an option in the FreqSum data format. + - Pileup Format now also records strand orientation 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.3.3 +version: 1.4.0 license: GPL-3 license-file: LICENSE maintainer: stephan.schiffels@mac.com diff --git a/src/SequenceFormats/FreqSum.hs b/src/SequenceFormats/FreqSum.hs --- a/src/SequenceFormats/FreqSum.hs +++ b/src/SequenceFormats/FreqSum.hs @@ -41,6 +41,7 @@ data FreqSumEntry = FreqSumEntry { fsChrom :: Chrom, -- ^The chromosome of the site fsPos :: Int, -- ^The position of the site + fsSnpId :: Maybe String, -- ^An optional parameter to take the snpId. This is not parsed from or printed to freqSum format but is used in internal conversions from Eigenstrat. fsRef :: Char, -- ^The reference allele fsAlt :: Char, -- ^The alternative allele fsCounts :: [Maybe Int] -- ^A list of allele counts in each group. Nothing denotes missing data. @@ -48,7 +49,7 @@ -- |This function converts a single freqSum entry to a printable freqSum line. freqSumEntryToText :: FreqSumEntry -> B.ByteString -freqSumEntryToText (FreqSumEntry chrom pos ref alt maybeCounts) = +freqSumEntryToText (FreqSumEntry chrom pos _ ref alt maybeCounts) = B.intercalate "\t" [B.pack (unChrom chrom), B.pack (show pos), B.singleton ref, B.singleton alt, countStr] <> "\n" where countStr = B.intercalate "\t" . map (B.pack . show . convertToNum) $ maybeCounts @@ -84,7 +85,7 @@ parseFreqSumEntry :: A.Parser FreqSumEntry parseFreqSumEntry = FreqSumEntry <$> (Chrom . B.unpack <$> A.takeTill isSpace) <* A.skipSpace <*> A.decimal <* - A.skipSpace <*> base <* A.skipSpace <*> baseOrDot <* A.skipSpace <*> counts <* A.endOfLine + A.skipSpace <*> pure Nothing <*> base <* A.skipSpace <*> baseOrDot <* A.skipSpace <*> counts <* A.endOfLine where counts = (parseMissing <|> parseCount) `A.sepBy` A.char '\t' parseMissing = A.string "-1" *> pure Nothing diff --git a/src/SequenceFormats/Pileup.hs b/src/SequenceFormats/Pileup.hs --- a/src/SequenceFormats/Pileup.hs +++ b/src/SequenceFormats/Pileup.hs @@ -1,5 +1,5 @@ {-# LANGUAGE OverloadedStrings #-} -module SequenceFormats.Pileup (readPileupFromStdIn, readPileupFromFile, PileupRow(..)) where +module SequenceFormats.Pileup (readPileupFromStdIn, readPileupFromFile, PileupRow(..), Strand(..)) where import Control.Monad.Catch (MonadThrow) import Control.Monad.IO.Class (MonadIO) @@ -12,10 +12,19 @@ import SequenceFormats.Utils (Chrom(..), word, readFileProd, consumeProducer) +-- |A datatype to represent the strand orientation of a single base. +data Strand = ForwardStrand | ReverseStrand deriving (Eq, Show) + -- |A datatype to represent a single pileup row for multiple individuals. -- The constructor arguments are: Chromosome, Position, Refererence Allelele, -- Pileup String per individual -data PileupRow = PileupRow Chrom Int Char [String] deriving (Eq, Show) +data PileupRow = PileupRow { + pileupChrom :: Chrom, -- ^The chromosome + pileupPos :: Int, -- ^The position + pileupRef :: Char, -- ^The reference base + pileupBases :: [String], -- ^The base string + pileupStrandInfo :: [[Strand]] + } deriving (Eq, Show) -- |Read a pileup-formatted file from StdIn, for reading from an -- external command `samtools mpileup`. @@ -36,10 +45,12 @@ -- for some reason, there is an M in the human reference at -- position 3:60830534 (both in hs37d5 and in hg19) _ <- A.space - entries <- parsePileupPerSample refA `A.sepBy1` + baseAndStrandEntries <- parsePileupPerSample refA `A.sepBy1` A.satisfy (\c -> c == ' ' || c == '\t') A.endOfLine - let ret = PileupRow (Chrom $ B.unpack chrom) pos refA entries + let baseStrings = map fst baseAndStrandEntries + strandInfoStrings = map snd baseAndStrandEntries + let ret = PileupRow (Chrom $ B.unpack chrom) pos refA baseStrings strandInfoStrings --trace (show ret) $ return ret return ret where @@ -47,13 +58,17 @@ processPileupEntry refA <$> A.decimal <* A.space <*> (B.unpack <$> word) <* A.space <* word -processPileupEntry :: Char -> Int -> String -> String +processPileupEntry :: Char -> Int -> String -> (String, [Strand]) processPileupEntry refA cov readBaseString = - if cov == 0 then "" else go readBaseString + if cov == 0 then ("", []) else + let res = go readBaseString + in (map fst res, map snd res) where go (x:xs) - | (x == '.' || x == ',') = refA : go xs - | x `elem` ("ACTGNactgn" :: String) = toUpper x : go xs + | x == '.' = (refA, ForwardStrand) : go xs + | 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 == '-') = diff --git a/src/SequenceFormats/Utils.hs b/src/SequenceFormats/Utils.hs --- a/src/SequenceFormats/Utils.hs +++ b/src/SequenceFormats/Utils.hs @@ -8,7 +8,7 @@ Chrom(..), word) where import Control.Error (readErr) -import Control.Exception (Exception) +import Control.Exception (Exception, throw) import Control.Monad.Catch (MonadThrow, throwM) import Control.Monad.Trans.Class (lift) import qualified Data.Attoparsec.ByteString.Char8 as A @@ -21,6 +21,12 @@ import qualified Pipes.Safe.Prelude as PS import System.IO (IOMode(..)) +-- |An exception type for parsing BioInformatic file formats. +data SeqFormatException = SeqFormatException String + deriving (Show, Eq) + +instance Exception SeqFormatException + -- |A wrapper datatype for Chromosome names. newtype Chrom = Chrom {unChrom :: String} deriving (Eq) @@ -31,20 +37,20 @@ -- |Ord instance for Chrom instance Ord Chrom where compare (Chrom c1) (Chrom c2) = - let c1' = if take 3 c1 == "chr" then drop 3 c1 else c1 - c2' = if take 3 c2 == "chr" then drop 3 c2 else c2 - in case (,) <$> readChrom c1' <*> readChrom c2' of - Left e -> error e + let [c1NoChr, c2NoChr] = map removeChr [c1, c2] + [c1XYMTconvert, c2XYMTconvert] = map convertXYMT [c1NoChr, c2NoChr] + in case (,) <$> readChrom c1XYMTconvert <*> readChrom c2XYMTconvert of + Left e -> throw $ SeqFormatException e Right (cn1, cn2) -> cn1 `compare` cn2 - -readChrom :: String -> Either String Int -readChrom c = readErr ("cannot parse chromosome " ++ c) $ c - --- |An exception type for parsing BioInformatic file formats. -data SeqFormatException = SeqFormatException String - deriving Show - -instance Exception SeqFormatException + where + removeChr c = if take 3 c == "chr" then drop 3 c else c + convertXYMT c = case c of + "X" -> "23" + "Y" -> "24" + "MT" -> "90" + n -> n + readChrom :: String -> Either String Int + readChrom c = readErr ("cannot parse chromosome " ++ c) $ c -- |A function to help with reporting parsing errors to stderr. Returns a clean Producer over the -- parsed datatype. diff --git a/src/SequenceFormats/VCF.hs b/src/SequenceFormats/VCF.hs --- a/src/SequenceFormats/VCF.hs +++ b/src/SequenceFormats/VCF.hs @@ -161,5 +161,5 @@ let alt = B.head . head . vcfAlt $ vcfEntry assertErr "Invalid Reference Allele" $ ref `elem` ['A', 'C', 'T', 'G', 'N'] assertErr "Invalid Alternative Allele" $ alt `elem` ['A', 'C', 'T', 'G', '.'] - return $ FreqSumEntry (vcfChrom vcfEntry) (vcfPos vcfEntry) ref alt dosages + return $ FreqSumEntry (vcfChrom vcfEntry) (vcfPos vcfEntry) (B.unpack <$> vcfId vcfEntry) ref alt dosages diff --git a/test/SequenceFormats/FreqSumSpec.hs b/test/SequenceFormats/FreqSumSpec.hs --- a/test/SequenceFormats/FreqSumSpec.hs +++ b/test/SequenceFormats/FreqSumSpec.hs @@ -49,11 +49,11 @@ mockDatFsEntries :: [FreqSumEntry] mockDatFsEntries = [ - FreqSumEntry (Chrom "11") 0 'A' 'C' [Just 1, Just 1, Just 1, Just 1, Just 1], - FreqSumEntry (Chrom "11") 100000 'A' 'G' [Just 2, Just 1, Just 0, Just 0, Just 0], - FreqSumEntry (Chrom "11") 200000 'A' 'T' [Just 0, Just 1, Just 1, Just 1, Just 1], - FreqSumEntry (Chrom "11") 300000 'C' 'A' [Just 2, Nothing, Just 1, Just 0, Just 0], - FreqSumEntry (Chrom "11") 400000 'G' 'A' [Just 0, Just 1, Just 1, Just 1, Just 1], - FreqSumEntry (Chrom "11") 500000 'T' 'A' [Just 2, Just 2, Just 1, Nothing, Just 1], - FreqSumEntry (Chrom "11") 600000 'G' 'T' [Just 0, Just 0, Just 1, Nothing, Nothing]] + FreqSumEntry (Chrom "11") 0 Nothing 'A' 'C' [Just 1, Just 1, Just 1, Just 1, Just 1], + FreqSumEntry (Chrom "11") 100000 Nothing 'A' 'G' [Just 2, Just 1, Just 0, Just 0, Just 0], + FreqSumEntry (Chrom "11") 200000 Nothing 'A' 'T' [Just 0, Just 1, Just 1, Just 1, Just 1], + FreqSumEntry (Chrom "11") 300000 Nothing 'C' 'A' [Just 2, Nothing, Just 1, Just 0, Just 0], + FreqSumEntry (Chrom "11") 400000 Nothing 'G' 'A' [Just 0, Just 1, Just 1, Just 1, Just 1], + FreqSumEntry (Chrom "11") 500000 Nothing 'T' 'A' [Just 2, Just 2, Just 1, Nothing, Just 1], + FreqSumEntry (Chrom "11") 600000 Nothing 'G' 'T' [Just 0, Just 0, Just 1, Nothing, Nothing]] diff --git a/test/SequenceFormats/PileupSpec.hs b/test/SequenceFormats/PileupSpec.hs --- a/test/SequenceFormats/PileupSpec.hs +++ b/test/SequenceFormats/PileupSpec.hs @@ -1,6 +1,6 @@ module SequenceFormats.PileupSpec (spec) where -import SequenceFormats.Pileup (readPileupFromFile, PileupRow(..)) +import SequenceFormats.Pileup (readPileupFromFile, PileupRow(..), Strand(..)) import SequenceFormats.Utils (Chrom(..)) import Control.Foldl (purely, list) @@ -19,6 +19,12 @@ mockDatPentries :: [PileupRow] mockDatPentries = [ - PileupRow (Chrom "1") 1000 'A' ["ACCT", "AAA", "ACCAACC"], - PileupRow (Chrom "1") 2000 'C' ["GGCA", "ACT", "ACCAACC"], - PileupRow (Chrom "2") 1000 'G' ["AGCT", "AAA", "ACCAACC"]] + PileupRow (Chrom "1") 1000 'A' ["AAACA", "AAAC", "AAAACCAACA"] + [[f, f, r, f, f], [f, f, r, r], [r, r, f, r, r, r, f, f, r, f]], + PileupRow (Chrom "1") 2000 'C' ["CCCA", "ACTCC", "CACACCCC"] + [[f, f, f, r], [f, f, r, f, f], [f, r, f, r, f, f, f, f]], + PileupRow (Chrom "2") 1000 'G' ["GGGGGGGCG", "GGG", "GGGGGG"] + [[r, r, f, r, r, f, f, f, f], [f, r, r], [f, r, r, f, r, r]]] + where + f = ForwardStrand + r = ReverseStrand diff --git a/test/SequenceFormats/UtilsSpec.hs b/test/SequenceFormats/UtilsSpec.hs --- a/test/SequenceFormats/UtilsSpec.hs +++ b/test/SequenceFormats/UtilsSpec.hs @@ -1,8 +1,9 @@ {-# LANGUAGE OverloadedStrings #-} module SequenceFormats.UtilsSpec (spec) where -import SequenceFormats.Utils (Chrom(..)) +import SequenceFormats.Utils (Chrom(..), SeqFormatException(..)) +import Control.Exception (evaluate) import Test.Hspec spec :: Spec @@ -18,3 +19,17 @@ Chrom "2" < Chrom "chr10" `shouldBe` True specify "chr2 should be smaller than chr10" $ Chrom "chr2" < Chrom "chr10" `shouldBe` True + specify "chr22 should be smaller than chrX" $ + Chrom "chr22" < Chrom "chrX" `shouldBe` True + specify "chrX should be smaller than chrY" $ + Chrom "chrX" < Chrom "chrY" `shouldBe` True + specify "chrY should be smaller than chrMT" $ + Chrom "chrY" < Chrom "chrMT" `shouldBe` True + specify "22 should be smaller than chrMT" $ + Chrom "22" < Chrom "chrMT" `shouldBe` True + specify "X should be smaller than chrMT" $ + Chrom "X" < Chrom "chrMT" `shouldBe` True + specify "chrSSS should throw" $ + evaluate (Chrom "chrSSS" < Chrom "chrMT") `shouldThrow` (==SeqFormatException "cannot parse chromosome SSS") + + diff --git a/test/SequenceFormats/VCFSpec.hs b/test/SequenceFormats/VCFSpec.hs --- a/test/SequenceFormats/VCFSpec.hs +++ b/test/SequenceFormats/VCFSpec.hs @@ -36,12 +36,12 @@ vcfRows !! 6 `shouldBe` vcf7 vcf1 :: VCFentry -vcf1 = VCFentry (Chrom "1") 10492 Nothing "C" ["T"] 15.0302 Nothing ["DP=28", "PV4=1,1,0.30985,1"] +vcf1 = VCFentry (Chrom "1") 10492 (Just "testId") "C" ["T"] 15.0302 Nothing ["DP=28", "PV4=1,1,0.30985,1"] ["GT", "PL"] [["0/0", "0,3,37"], ["0/0", "0,6,67"], ["0/1", "51,0,28"], ["0/0", "0,54,255"], ["0/0", "0,9,83"]] vcf7 :: VCFentry -vcf7 = VCFentry (Chrom "2") 30923 (Just "rs12345") "G" [] 110.112 Nothing ["DP=5", "FQ=-28.9619"] +vcf7 = VCFentry (Chrom "2") 30923 Nothing "G" [] 110.112 Nothing ["DP=5", "FQ=-28.9619"] ["GT", "PL"] [["1/1", "0,0,0"], ["1/1", "0,0,0"], ["1/1", "40,6,0"], ["1/1", "105,9,0"], ["1/1", "0,0,0"]] testGetGenotypes :: Spec @@ -74,7 +74,7 @@ testVcfToFreqsumEntry :: Spec testVcfToFreqsumEntry = describe "vcfToFreqsumEntry" $ it "should convert correctly" $ do - let r = Right (FreqSumEntry (Chrom "1") 10492 'C' 'T' [Just 0, Just 0, Just 1, Just 0, Just 0]) + let r = Right (FreqSumEntry (Chrom "1") 10492 (Just "testId") 'C' 'T' [Just 0, Just 0, Just 1, Just 0, Just 0]) vcfToFreqSumEntry vcf1 `shouldBe` r testIsBiallelicSnp :: Spec diff --git a/testDat/example.vcf b/testDat/example.vcf --- a/testDat/example.vcf +++ b/testDat/example.vcf @@ -18,10 +18,10 @@ ##bcftools_callVersion=1.3+htslib-1.3 ##bcftools_callCommand=call -c -v #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT 12880A 12881A 12883A 12884A 12885A -1 10492 . C T 15.0302 . DP=28;PV4=1,1,0.30985,1 GT:PL 0/0:0,3,37 0/0:0,6,67 0/1:51,0,28 0/0:0,54,255 0/0:0,9,83 +1 10492 testId C T 15.0302 . DP=28;PV4=1,1,0.30985,1 GT:PL 0/0:0,3,37 0/0:0,6,67 0/1:51,0,28 0/0:0,54,255 0/0:0,9,83 1 14907 . A G 146.092 . DP=55.526787;AC1=5;DP4=19,16,8,12;PV4=0.403198,0.000983733,1,0.0348559 GT:PL 0/1:0,0,0 0/0:0,9,42 0/1:51,0,0 0/1:70,0,123 1/1:68,5,0 1 14930 . A G 999 . DP=61.545999;AC1=6;;FQ=999;PV4=0.602167,0.240103,1,1 GT:PL 0/1:0,3,20 0/1:22,0,47 0/1:70,0,26 0/1:188,0,121 1/1:80,12,0 2 14933 . G A 51.2451 . DP=60.215796;PV4=0.49427,1,1,1 GT:PL 0/0:0,6,39 0/1:25,0,48 0/0:0,18,118 0/1:69,0,163 0/0:0,12,80 2 16495 . G C 7.29715 . DP=47.27107;AC1=2;PV4=0.540781,0.244272,1,0.0352549 GT:PL 0/0:0,9,9 0/1:23,3,0 0/0:0,4,9 0/1:27,0,46 0/0:0,9,21 2 20144 . G A 15.7039 . DP=121;AC1=3;PV4=1,1,1,1 GT:PL 0/0:0,5,81 0/1:33,0,47 0/1:9,0,21 0/0:0,182,122 0/1:18,0,27 -2 30923 rs12345 G . 110.112 . DP=5;FQ=-28.9619 GT:PL 1/1:0,0,0 1/1:0,0,0 1/1:40,6,0 1/1:105,9,0 1/1:0,0,0 +2 30923 . G . 110.112 . DP=5;FQ=-28.9619 GT:PL 1/1:0,0,0 1/1:0,0,0 1/1:40,6,0 1/1:105,9,0 1/1:0,0,0