packages feed

sequence-formats 1.3.3 → 1.4.0

raw patch · 11 files changed

+92/−44 lines, 11 files

Files

Changelog.md view
@@ -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
sequence-formats.cabal view
@@ -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
src/SequenceFormats/FreqSum.hs view
@@ -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
src/SequenceFormats/Pileup.hs view
@@ -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 == '-') =
src/SequenceFormats/Utils.hs view
@@ -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.
src/SequenceFormats/VCF.hs view
@@ -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
test/SequenceFormats/FreqSumSpec.hs view
@@ -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]]
test/SequenceFormats/PileupSpec.hs view
@@ -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
test/SequenceFormats/UtilsSpec.hs view
@@ -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")+    +
test/SequenceFormats/VCFSpec.hs view
@@ -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
testDat/example.vcf view
@@ -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