sequence-formats 1.3.2.1 → 1.3.3
raw patch · 21 files changed
+423/−325 lines, 21 filesdep +hspec
Dependencies added: hspec
Files
- Changelog.md +2/−0
- sequence-formats.cabal +13/−10
- src-tests/SequenceFormats/Eigenstrat/Test.hs +0/−79
- src-tests/SequenceFormats/Fasta/Test.hs +0/−17
- src-tests/SequenceFormats/FreqSum/Test.hs +0/−49
- src-tests/SequenceFormats/RareAlleleHistogram/Test.hs +0/−40
- src-tests/SequenceFormats/Utils/Test.hs +0/−12
- src-tests/SequenceFormats/VCF/Test.hs +0/−61
- src-tests/tests.hs +0/−48
- src/SequenceFormats/Eigenstrat.hs +2/−5
- src/SequenceFormats/Pileup.hs +62/−0
- src/SequenceFormats/Utils.hs +5/−2
- src/SequenceFormats/VCF.hs +2/−2
- test/SequenceFormats/EigenstratSpec.hs +81/−0
- test/SequenceFormats/FastaSpec.hs +23/−0
- test/SequenceFormats/FreqSumSpec.hs +59/−0
- test/SequenceFormats/PileupSpec.hs +24/−0
- test/SequenceFormats/RareAlleleHistogramSpec.hs +44/−0
- test/SequenceFormats/UtilsSpec.hs +20/−0
- test/SequenceFormats/VCFSpec.hs +85/−0
- test/Spec.hs +1/−0
Changelog.md view
@@ -25,3 +25,5 @@ V 1.3.2: Added testDat to Cabal file to make tests work off the tarball. 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.
sequence-formats.cabal view
@@ -1,6 +1,6 @@ cabal-version: >=1.10 name: sequence-formats-version: 1.3.2.1+version: 1.3.3 license: GPL-3 license-file: LICENSE maintainer: stephan.schiffels@mac.com@@ -30,6 +30,7 @@ SequenceFormats.VCF SequenceFormats.Eigenstrat SequenceFormats.Utils+ SequenceFormats.Pileup hs-source-dirs: src default-language: Haskell2010 build-depends:@@ -50,15 +51,16 @@ test-suite sequenceFormatTests type: exitcode-stdio-1.0- main-is: tests.hs- hs-source-dirs: src-tests+ main-is: Spec.hs+ hs-source-dirs: test other-modules:- SequenceFormats.Eigenstrat.Test- SequenceFormats.Fasta.Test- SequenceFormats.FreqSum.Test- SequenceFormats.RareAlleleHistogram.Test- SequenceFormats.Utils.Test- SequenceFormats.VCF.Test+ SequenceFormats.EigenstratSpec+ SequenceFormats.FastaSpec+ SequenceFormats.FreqSumSpec+ SequenceFormats.RareAlleleHistogramSpec+ SequenceFormats.UtilsSpec+ SequenceFormats.VCFSpec+ SequenceFormats.PileupSpec default-language: Haskell2010 build-depends: base >=4.12.0.0,@@ -71,4 +73,5 @@ transformers >=0.5.6.2, tasty-hunit >=0.10.0.2, bytestring >=0.10.8.2,- containers >=0.6.0.1+ containers >=0.6.0.1,+ hspec >=2.7.1
− src-tests/SequenceFormats/Eigenstrat/Test.hs
@@ -1,79 +0,0 @@-{-# LANGUAGE OverloadedStrings #-}-module SequenceFormats.Eigenstrat.Test (bimReadTest, eigenstratReadTest, eigenstratWriteTest) where--import Control.Foldl (purely, list)-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 SequenceFormats.Eigenstrat (readEigenstrat, writeEigenstrat, readBimFile,- EigenstratSnpEntry(..), EigenstratIndEntry(..), GenoLine, Sex(..), GenoEntry(..))-import SequenceFormats.Utils (Chrom(..))-import Test.Tasty.HUnit (Assertion, assertEqual)--bimReadTest :: Assertion-bimReadTest = do- let esSnpProd = readBimFile "testDat/example.bim"- esSnpDat <- runSafeT $ purely P.fold list esSnpProd- assertEqual "" testDatEigenstratSnp esSnpDat --eigenstratReadTest :: Assertion-eigenstratReadTest = do- let esSnpFile = "testDat/example.snp"- esIndFile = "testDat/example.ind"- esGenoFile = "testDat/example.eigenstratgeno"- (indEntries, esProd) <- runSafeT $ readEigenstrat esGenoFile esSnpFile esIndFile- assertEqual "eigenstratReadTest_assertIndEntries" testDatEigenstratInd indEntries- snpGenoEntries <- runSafeT $ purely P.fold list esProd- assertEqual "eigenstratReadTest_assertSnpEntries" testDatEigenstratSnp (map fst snpGenoEntries)- assertEqual "eigenstratReadTest_assertGenoEntries" testDatEigenstratGeno- (map snd snpGenoEntries)--eigenstratWriteTest :: Assertion-eigenstratWriteTest = do- let tmpGeno = "/tmp/eigenstratWriteTest.geno"- tmpSnp = "/tmp/eigenstratWriteTest.snp"- tmpInd = "/tmp/eigenstratWriteTest.ind"- testDatSnpProd = each testDatEigenstratSnp- testDatGenoProd = each testDatEigenstratGeno- testDatJointProd = P.zip testDatSnpProd testDatGenoProd- liftIO . runSafeT . runEffect $- testDatJointProd >-> writeEigenstrat tmpGeno tmpSnp tmpInd testDatEigenstratInd- (indEntries, esProd) <- liftIO . runSafeT $ readEigenstrat tmpGeno tmpSnp tmpInd- liftIO $ assertEqual "eigenstratWriteTest_assertIndEntries" testDatEigenstratInd indEntries- snpGenoEntries <- liftIO . runSafeT $ purely P.fold list esProd- liftIO $ assertEqual "eigenstratWriteTest_assertIndEntries" testDatEigenstratSnp- (map fst snpGenoEntries)- liftIO $ assertEqual "eigenstratWriteTest_assertIndEntries" testDatEigenstratGeno- (map snd snpGenoEntries)--testDatEigenstratSnp :: [EigenstratSnpEntry]-testDatEigenstratSnp = [- EigenstratSnpEntry (Chrom "11") 0 0.000000 "rs0000" 'A' 'C',- EigenstratSnpEntry (Chrom "11") 100000 0.001000 "rs1111" 'A' 'G',- EigenstratSnpEntry (Chrom "11") 200000 0.002000 "rs2222" 'A' 'T',- EigenstratSnpEntry (Chrom "11") 300000 0.003000 "rs3333" 'C' 'A',- EigenstratSnpEntry (Chrom "11") 400000 0.004000 "rs4444" 'G' 'A',- EigenstratSnpEntry (Chrom "11") 500000 0.005000 "rs5555" 'T' 'A',- EigenstratSnpEntry (Chrom "11") 600000 0.006000 "rs6666" 'G' 'T']---testDatEigenstratInd :: [EigenstratIndEntry]-testDatEigenstratInd = [- EigenstratIndEntry "SAMPLE0" Female "Case",- EigenstratIndEntry "SAMPLE1" Male "Case",- EigenstratIndEntry "SAMPLE2" Female "Control",- EigenstratIndEntry "SAMPLE3" Male "Control",- EigenstratIndEntry "SAMPLE4" Female "Control"]--testDatEigenstratGeno :: [GenoLine]-testDatEigenstratGeno = [- fromList [Het, Het, Het, HomAlt, HomAlt],- fromList [HomAlt, Het, HomRef, Het, HomRef],- fromList [HomRef, Het, Het, HomAlt, Het],- fromList [HomAlt, Missing, Het, HomRef, HomRef],- fromList [HomRef, Het, Het, HomAlt, HomAlt],- fromList [HomAlt, HomAlt, Het, Missing, Het],- fromList [HomRef, HomRef, Het, Missing, Missing]]-
− src-tests/SequenceFormats/Fasta/Test.hs
@@ -1,17 +0,0 @@-{-# LANGUAGE OverloadedStrings #-}-module SequenceFormats.Fasta.Test (testFastaRead) where--import qualified Data.ByteString.Char8 as BS-import qualified Pipes.Prelude as P-import SequenceFormats.Fasta (loadFastaChrom)-import SequenceFormats.Utils (Chrom(..))-import System.IO (withFile, IOMode(..))-import Test.Tasty.HUnit (Assertion, assertEqual)--testFastaRead :: Assertion-testFastaRead = withFile "testDat/example.fasta" ReadMode $ \h -> do- fastaProd <- loadFastaChrom h (Chrom "chr3")- fastaString <- BS.unpack . BS.concat <$> P.toListM fastaProd- assertEqual "fastaReadTest_assertFastaEntries" testString fastaString- where- testString = "ACGACGACGACGGGGTTTAAAAAGGGTTTCCTCTCTCTCTGGG"
− src-tests/SequenceFormats/FreqSum/Test.hs
@@ -1,49 +0,0 @@-{-# LANGUAGE OverloadedStrings #-}-module SequenceFormats.FreqSum.Test (fsReadTest, fsWriteTest) where--import SequenceFormats.FreqSum (readFreqSumFile, printFreqSumFile, FreqSumEntry(..), - FreqSumHeader(..))-import SequenceFormats.Utils (Chrom(..))--import Control.Foldl (purely, list)-import Control.Monad.IO.Class (liftIO)-import Pipes (each, runEffect, (>->))-import qualified Pipes.Prelude as P-import Pipes.Safe (runSafeT)-import Test.Tasty.HUnit (Assertion, assertEqual)--fsReadTest :: Assertion-fsReadTest = runSafeT $ do- let fsFile = "testDat/example.freqsum"- (fsHeader, fsProd) <- readFreqSumFile fsFile- liftIO $ assertEqual "fsReadTest_assertIndEntries" testDatFsHeader fsHeader- fsEntries <- purely P.fold list fsProd- liftIO $ assertEqual "fsReadTest_assertFsEntries" testDatFsEntries fsEntries--fsWriteTest :: Assertion-fsWriteTest = do- let fn = "/tmp/freqSumWriteTest.txt"- testDatFsProd = each testDatFsEntries- runSafeT . runEffect $ testDatFsProd >-> printFreqSumFile fn testDatFsHeader- runSafeT $ do- (fsHeader, fsProd) <- readFreqSumFile fn- liftIO $ assertEqual "fsWriteTest_assertIndEntries" testDatFsHeader fsHeader- fsEntries <- purely P.fold list fsProd- liftIO $ assertEqual "fsWriteTest_assertFsEntries" testDatFsEntries fsEntries--testDatFsHeader :: FreqSumHeader -testDatFsHeader = FreqSumHeader names numbers- where- names = ["SAMPLE0", "SAMPLE1", "SAMPLE2", "SAMPLE3", "SAMPLE4"]- numbers = [2, 2, 2, 1, 1]--testDatFsEntries :: [FreqSumEntry]-testDatFsEntries = [- 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]]-
− src-tests/SequenceFormats/RareAlleleHistogram/Test.hs
@@ -1,40 +0,0 @@-{-# LANGUAGE OverloadedStrings #-}-module SequenceFormats.RareAlleleHistogram.Test (testReadHistogram, testWriteHistogram) where--import SequenceFormats.RareAlleleHistogram (RareAlleleHistogram(..), readHistogram, writeHistogramFile)--import Control.Monad.IO.Class (liftIO)-import qualified Data.Map as Map-import Test.Tasty.HUnit (Assertion, assertEqual)--testReadHistogram :: Assertion-testReadHistogram = do- hist <- readHistogram "testDat/example.histogram.txt"- assertEqual "readHistogramTest" hist testHistogramDat--testWriteHistogram :: Assertion-testWriteHistogram = do- let fn = "/tmp/histogramWriteTest.txt"- writeHistogramFile fn testHistogramDat- hist <- liftIO $ readHistogram fn- liftIO $ assertEqual "writeHistogramTest" hist testHistogramDat--testHistogramDat :: RareAlleleHistogram-testHistogramDat = RareAlleleHistogram names nVec 1 10 [] [] 1146826657 counts jnEstimates- where- names = ["EUR", "SEA", "SIB", "CHK", "SAM"]- nVec = [66, 42, 44, 8, 28]- counts = Map.fromList [- ([1,0,0,0,0], 773148),- ([0,1,0,0,0], 527207),- ([0,0,1,0,0], 368640),- ([0,0,0,0,1], 213918),- ([2,0,0,0,0], 158795),- ([0,2,0,0,0], 73241)]- jnEstimates = Just $ Map.fromList [- ([1,0,0,0,0], (6.7e-4, 3.6e-6)),- ([0,1,0,0,0], (4.5e-4, 2.5e-6)),- ([0,0,1,0,0], (3.2e-4, 2.4e-6)),- ([0,0,0,0,1], (1.8e-4, 3.6e-6)),- ([2,0,0,0,0], (1.3e-4, 1.4e-6)),- ([0,2,0,0,0], (6.3e-5, 7.3e-7))]
− src-tests/SequenceFormats/Utils/Test.hs
@@ -1,12 +0,0 @@-{-# LANGUAGE OverloadedStrings #-}-module SequenceFormats.Utils.Test (testChromOrder) where--import SequenceFormats.Utils (Chrom(..))-import Test.Tasty.HUnit (Assertion, assertBool)--testChromOrder :: Assertion-testChromOrder = do- assertBool "testChromSmaller" (Chrom "2" < Chrom "10")- assertBool "testChromSmaller2" (Chrom "chr2" < Chrom "10")- assertBool "testChromSmaller3" (Chrom "2" < Chrom "chr10")- assertBool "testChromSmaller3" (Chrom "chr2" < Chrom "chr10")
− src-tests/SequenceFormats/VCF/Test.hs
@@ -1,61 +0,0 @@-{-# LANGUAGE OverloadedStrings #-}-module SequenceFormats.VCF.Test (testReadVCF, testGetGenotypes, testGetDosages, testIsTransversionSnp, testIsBiallelicSnp, testVcfToFreqsumEntry) where--import Control.Foldl (list, purely)-import Control.Monad.IO.Class (liftIO)-import Pipes.Prelude (fold)-import Pipes.Safe (runSafeT)-import SequenceFormats.FreqSum (FreqSumEntry(..))-import SequenceFormats.Utils (Chrom(..))-import SequenceFormats.VCF (readVCFfromFile, getGenotypes, getDosages, isTransversionSnp, vcfToFreqSumEntry, isBiallelicSnp, VCFheader(..), VCFentry(..))-import Test.Tasty.HUnit (Assertion, assertEqual)--testReadVCF :: Assertion-testReadVCF = runSafeT $ do- (vcfH, vcfProd) <- readVCFfromFile "testDat/example.vcf"- let vcfHc = vcfHeaderComments vcfH- liftIO $ assertEqual "headerLine1" "##fileformat=VCFv4.2" (vcfHc !! 0)- liftIO $ assertEqual "headerLine19" "##bcftools_callCommand=call -c -v" (vcfHc !! 18)- let n = vcfSampleNames vcfH- liftIO $ assertEqual "sample names" ["12880A", "12881A", "12883A", "12884A", "12885A"] n- rows <- purely fold list vcfProd- liftIO $ assertEqual "vcf1" vcf1 (rows !! 0)- liftIO $ assertEqual "vcf7" vcf7 (rows !! 6)--vcf1 :: VCFentry-vcf1 = VCFentry (Chrom "1") 10492 Nothing "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"]- ["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 :: Assertion-testGetGenotypes = do- assertEqual "getGenotypesRight" (Right ["0/0", "0/0", "0/1", "0/0", "0/0"]) (getGenotypes vcf1)- assertEqual "getGenotypesLeft" (Left "GT format field not found") (getGenotypes (vcf1 {vcfFormatString=["PL"]}))--testGetDosages :: Assertion-testGetDosages = do- assertEqual "getDosages1" (Right [Just 0, Just 0, Just 1, Just 0, Just 0]) (getDosages vcf1)- assertEqual "getDosages2" (Right [Just 0, Just 0, Nothing, Just 0, Just 0])- (getDosages (vcf1 {vcfGenotypeInfo=[["0/0", "0,3,37"], ["0/0", "0,6,67"], [".", "51,0,28"], ["0/0", "0,54,255"],- ["0/0", "0,9,83"]]}))--testIsTransversionSnp :: Assertion-testIsTransversionSnp = do- assertEqual "testTriAllelic" False (isTransversionSnp "A" ["C", "T"])- assertEqual "testTransition1" False (isTransversionSnp "A" ["G"])- assertEqual "testTransition2" False (isTransversionSnp "T" ["C"])- assertEqual "testTransversion1" True (isTransversionSnp "T" ["A"])- assertEqual "testTransversion2" True (isTransversionSnp "C" ["G"])--testVcfToFreqsumEntry :: Assertion-testVcfToFreqsumEntry =- assertEqual "vcfToFreqsumValid" (Right (FreqSumEntry (Chrom "1") 10492 'C' 'T' [Just 0, Just 0, Just 1, Just 0, Just 0])) (vcfToFreqSumEntry vcf1)--testIsBiallelicSnp :: Assertion-testIsBiallelicSnp = do- assertEqual "testTriAllelic" False (isBiallelicSnp "A" ["C", "T"])- assertEqual "testBiallelic" True (isBiallelicSnp "A" ["T"])
− src-tests/tests.hs
@@ -1,48 +0,0 @@-{-# LANGUAGE OverloadedStrings #-}--import SequenceFormats.RareAlleleHistogram.Test (testWriteHistogram, testReadHistogram)-import SequenceFormats.Eigenstrat.Test (eigenstratReadTest, bimReadTest, eigenstratWriteTest)-import SequenceFormats.FreqSum.Test (fsReadTest, fsWriteTest)-import SequenceFormats.Fasta.Test (testFastaRead)-import SequenceFormats.Utils.Test (testChromOrder)-import SequenceFormats.VCF.Test (testReadVCF, testGetGenotypes, testGetDosages, testIsTransversionSnp, testIsBiallelicSnp, testVcfToFreqsumEntry)-import Test.Tasty (TestTree, testGroup, defaultMain)-import Test.Tasty.HUnit (testCase)--main :: IO ()-main = defaultMain allTests- where- allTests = testGroup "Tests" [eigenstratTests, freqSumTests, fastaTests, rahTests, chromTests,- vcfTests]--eigenstratTests :: TestTree-eigenstratTests = testGroup "Eigenstrat Tests" [- testCase "Eigenstrat Read Test" eigenstratReadTest,- testCase "Bim Read Test" bimReadTest,- testCase "Eigenstrat Write Test" eigenstratWriteTest]--freqSumTests :: TestTree-freqSumTests = testGroup "FreqSum Tests" [- testCase "FreqSum Read Test" fsReadTest,- testCase "FreqSum Write Test" fsWriteTest]--fastaTests :: TestTree-fastaTests = testGroup "Fasta Tests" [testCase "Fasta Read Test" testFastaRead]--rahTests :: TestTree-rahTests = testGroup "Histogram Tests" [- testCase "Histogram Read Test" testReadHistogram,- testCase "Histogram Write Test" testWriteHistogram]--chromTests :: TestTree-chromTests = testGroup "Chrom tests" [- testCase "Chromosome ordering tests" testChromOrder]--vcfTests :: TestTree-vcfTests = testGroup "VCF tests" [- testCase "VCF read test" testReadVCF,- testCase "VCF getGenotypes" testGetGenotypes,- testCase "VCF getDosages" testGetDosages,- testCase "VCF isTransitionSnp" testIsTransversionSnp,- testCase "VCF isBiallelicSnp" testIsBiallelicSnp,- testCase "VCF testVcfToFreqsumEntry" testVcfToFreqsumEntry]
src/SequenceFormats/Eigenstrat.hs view
@@ -10,7 +10,8 @@ readEigenstrat, writeEigenstrat, writeEigenstratIndFile, writeEigenstratSnp, writeBim, writeEigenstratGeno) where -import SequenceFormats.Utils (consumeProducer, SeqFormatException(..), Chrom(..), readFileProd)+import SequenceFormats.Utils (consumeProducer, SeqFormatException(..),+ Chrom(..), readFileProd, word) import Control.Applicative ((<|>)) import Control.Exception (throw)@@ -18,7 +19,6 @@ import Control.Monad.Catch (MonadThrow) import Control.Monad.IO.Class (MonadIO, liftIO) import qualified Data.Attoparsec.ByteString.Char8 as A-import Data.Char (isSpace) import Data.Vector (Vector, fromList, toList) import qualified Data.ByteString.Char8 as B import Pipes (Producer, Pipe, (>->), for, cat, yield, Consumer)@@ -74,9 +74,6 @@ void A.endOfLine return $ EigenstratSnpEntry (Chrom (B.unpack chrom)) pos geneticPos snpId_ ref alt -word :: A.Parser B.ByteString-word = A.takeTill isSpace- eigenstratIndParser :: A.Parser EigenstratIndEntry eigenstratIndParser = do A.skipMany A.space
+ src/SequenceFormats/Pileup.hs view
@@ -0,0 +1,62 @@+{-# LANGUAGE OverloadedStrings #-}+module SequenceFormats.Pileup (readPileupFromStdIn, readPileupFromFile, PileupRow(..)) where++import Control.Monad.Catch (MonadThrow)+import Control.Monad.IO.Class (MonadIO)+import qualified Data.Attoparsec.ByteString.Char8 as A+import qualified Data.ByteString.Char8 as B+import Data.Char (toUpper)+import Pipes (Producer)+import qualified Pipes.ByteString as PB+import Pipes.Safe (MonadSafe)++import SequenceFormats.Utils (Chrom(..), word, readFileProd, consumeProducer)++-- |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)++-- |Read a pileup-formatted file from StdIn, for reading from an+-- external command `samtools mpileup`.+readPileupFromStdIn :: (MonadIO m, MonadThrow m) => Producer PileupRow m ()+readPileupFromStdIn = consumeProducer pileupParser PB.stdin++-- |Read pileup from a file.+readPileupFromFile :: (MonadSafe m) => FilePath -> Producer PileupRow m ()+readPileupFromFile = consumeProducer pileupParser . readFileProd++pileupParser :: A.Parser PileupRow+pileupParser = do+ chrom <- word+ _ <- A.space+ pos <- A.decimal+ _ <- A.space+ refA <- A.satisfy (A.inClass "ACTGNactgnM")+ -- 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`+ A.satisfy (\c -> c == ' ' || c == '\t')+ A.endOfLine+ let ret = PileupRow (Chrom $ B.unpack chrom) pos refA entries+ --trace (show ret) $ return ret+ return ret+ where+ parsePileupPerSample refA =+ processPileupEntry refA <$> A.decimal <* A.space <*> (B.unpack <$> word) <*+ A.space <* word++processPileupEntry :: Char -> Int -> String -> String+processPileupEntry refA cov readBaseString =+ if cov == 0 then "" else go readBaseString+ where+ go (x:xs)+ | (x == '.' || x == ',') = refA : go xs+ | x `elem` ("ACTGNactgn" :: String) = toUpper x : go xs+ | x `elem` ("$*" :: String) = go xs+ | x == '^' = go (drop 1 xs)+ | (x == '+' || x == '-') =+ let [(num, rest)] = reads xs in go (drop num rest)+ | otherwise = error $ "cannot parse read base string: " ++ (x:xs)+ go [] = []
src/SequenceFormats/Utils.hs view
@@ -5,14 +5,15 @@ module SequenceFormats.Utils (liftParsingErrors, consumeProducer, readFileProd, SeqFormatException(..),- Chrom(..)) where+ Chrom(..), word) where import Control.Error (readErr) import Control.Exception (Exception) import Control.Monad.Catch (MonadThrow, throwM) import Control.Monad.Trans.Class (lift)-import qualified Data.ByteString.Char8 as B import qualified Data.Attoparsec.ByteString.Char8 as A+import qualified Data.ByteString.Char8 as B+import Data.Char (isSpace) import Pipes (Producer, next) import Pipes.Attoparsec (ParsingError(..), parsed) import qualified Pipes.ByteString as PB@@ -66,3 +67,5 @@ readFileProd :: (PS.MonadSafe m) => FilePath -> Producer B.ByteString m () readFileProd f = PS.withFile f ReadMode PB.fromHandle +word :: A.Parser B.ByteString+word = A.takeTill isSpace
src/SequenceFormats/VCF.hs view
@@ -15,7 +15,8 @@ vcfToFreqSumEntry, isBiallelicSnp) where -import SequenceFormats.Utils (consumeProducer, Chrom(..), readFileProd, SeqFormatException(..))+import SequenceFormats.Utils (consumeProducer, Chrom(..),+ readFileProd, SeqFormatException(..), word) import SequenceFormats.FreqSum (FreqSumEntry(..)) import Control.Applicative ((<|>))@@ -97,7 +98,6 @@ vcfEntryParserTruncated = VCFentry <$> (Chrom . B.unpack <$> word) <* sp <*> A.decimal <* sp <*> parseId <* sp <*> word <* sp <*> parseAlternativeAlleles <* sp <*> A.double <* sp <*> parseFilter <* sp <*> parseInfoFields <*> pure [] <*> pure [] <* A.endOfLine- word = A.takeTill isSpace sp = A.satisfy (\c -> c == ' ' || c == '\t') parseId = (parseDot *> pure Nothing) <|> (Just <$> word) parseDot = A.char '.'
+ test/SequenceFormats/EigenstratSpec.hs view
@@ -0,0 +1,81 @@+{-# LANGUAGE OverloadedStrings #-}+module SequenceFormats.EigenstratSpec (spec) where++import Control.Foldl (purely, list)+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 SequenceFormats.Eigenstrat (readEigenstrat, writeEigenstrat, readBimFile,+ EigenstratSnpEntry(..), EigenstratIndEntry(..), GenoLine, Sex(..), GenoEntry(..))+import SequenceFormats.Utils (Chrom(..))+import Test.Hspec++spec :: Spec+spec = do+ testReadBimFile+ testReadEigenstrat+ testWriteEigenstrat++mockDatEigenstratSnp :: [EigenstratSnpEntry]+mockDatEigenstratSnp = [+ EigenstratSnpEntry (Chrom "11") 0 0.000000 "rs0000" 'A' 'C',+ EigenstratSnpEntry (Chrom "11") 100000 0.001000 "rs1111" 'A' 'G',+ EigenstratSnpEntry (Chrom "11") 200000 0.002000 "rs2222" 'A' 'T',+ EigenstratSnpEntry (Chrom "11") 300000 0.003000 "rs3333" 'C' 'A',+ EigenstratSnpEntry (Chrom "11") 400000 0.004000 "rs4444" 'G' 'A',+ EigenstratSnpEntry (Chrom "11") 500000 0.005000 "rs5555" 'T' 'A',+ EigenstratSnpEntry (Chrom "11") 600000 0.006000 "rs6666" 'G' 'T']++testReadBimFile :: Spec+testReadBimFile = describe "readBimFile" $+ it "should read a BIM file correctly" $ do+ let esSnpProd = readBimFile "testDat/example.bim"+ (runSafeT $ purely P.fold list esSnpProd) `shouldReturn` mockDatEigenstratSnp++mockDatEigenstratInd :: [EigenstratIndEntry]+mockDatEigenstratInd = [+ EigenstratIndEntry "SAMPLE0" Female "Case",+ EigenstratIndEntry "SAMPLE1" Male "Case",+ EigenstratIndEntry "SAMPLE2" Female "Control",+ EigenstratIndEntry "SAMPLE3" Male "Control",+ EigenstratIndEntry "SAMPLE4" Female "Control"]++mockDatEigenstratGeno :: [GenoLine]+mockDatEigenstratGeno = [+ fromList [Het, Het, Het, HomAlt, HomAlt],+ fromList [HomAlt, Het, HomRef, Het, HomRef],+ fromList [HomRef, Het, Het, HomAlt, Het],+ fromList [HomAlt, Missing, Het, HomRef, HomRef],+ fromList [HomRef, Het, Het, HomAlt, HomAlt],+ fromList [HomAlt, HomAlt, Het, Missing, Het],+ fromList [HomRef, HomRef, Het, Missing, Missing]]+ +testReadEigenstrat :: Spec+testReadEigenstrat = describe "readEigenstrat" $ do+ it "should read the correct eigenstrat file" $ do+ let esSnpFile = "testDat/example.snp"+ esIndFile = "testDat/example.ind"+ esGenoFile = "testDat/example.eigenstratgeno"+ (indEntries, esProd) <- runSafeT $ readEigenstrat esGenoFile esSnpFile esIndFile+ indEntries `shouldBe` mockDatEigenstratInd + snpGenoEntries <- runSafeT $ purely P.fold list esProd+ (map fst snpGenoEntries) `shouldBe` mockDatEigenstratSnp + (map snd snpGenoEntries) `shouldBe` mockDatEigenstratGeno++testWriteEigenstrat :: Spec+testWriteEigenstrat = describe "writeEigenstrat" $ do+ it "should write and read back eigenstrat data correctly" $ do+ let tmpGeno = "/tmp/eigenstratWriteTest.geno"+ tmpSnp = "/tmp/eigenstratWriteTest.snp"+ tmpInd = "/tmp/eigenstratWriteTest.ind"+ testDatSnpProd = each mockDatEigenstratSnp+ testDatGenoProd = each mockDatEigenstratGeno+ testDatJointProd = P.zip testDatSnpProd testDatGenoProd+ liftIO . runSafeT . runEffect $+ testDatJointProd >-> writeEigenstrat tmpGeno tmpSnp tmpInd mockDatEigenstratInd+ (indEntries, esProd) <- liftIO . runSafeT $ readEigenstrat tmpGeno tmpSnp tmpInd+ indEntries `shouldBe` mockDatEigenstratInd + snpGenoEntries <- liftIO . runSafeT $ purely P.fold list esProd+ (map fst snpGenoEntries) `shouldBe` mockDatEigenstratSnp
+ test/SequenceFormats/FastaSpec.hs view
@@ -0,0 +1,23 @@+{-# LANGUAGE OverloadedStrings #-}+module SequenceFormats.FastaSpec (spec) where++import qualified Data.ByteString.Char8 as BS+import qualified Pipes.Prelude as P+import SequenceFormats.Fasta (loadFastaChrom)+import SequenceFormats.Utils (Chrom(..))+import System.IO (withFile, IOMode(..))+import Test.Hspec++spec :: Spec+spec = do+ testLoadFastaChrom++testLoadFastaChrom :: Spec+testLoadFastaChrom = around (withFile "testDat/example.fasta" ReadMode) $+ describe "loadFastaChrom" $+ it "should read the correct fasta string from file" $ \h -> do+ fastaProd <- loadFastaChrom h (Chrom "chr3")+ fastaString <- BS.unpack . BS.concat <$> P.toListM fastaProd+ fastaString `shouldBe` testString+ where+ testString = "ACGACGACGACGGGGTTTAAAAAGGGTTTCCTCTCTCTCTGGG"
+ test/SequenceFormats/FreqSumSpec.hs view
@@ -0,0 +1,59 @@+{-# LANGUAGE OverloadedStrings #-}+module SequenceFormats.FreqSumSpec (spec) where++import SequenceFormats.FreqSum (readFreqSumFile, printFreqSumFile, FreqSumEntry(..), + FreqSumHeader(..))+import SequenceFormats.Utils (Chrom(..))++import Control.Foldl (purely, list)+import Pipes (each, runEffect, (>->))+import qualified Pipes.Prelude as P+import Pipes.Safe (runSafeT)+import Test.Hspec++spec :: Spec+spec = do+ testReadFreqSumFile+ testPrintFreqSumFile++testReadFreqSumFile :: Spec+testReadFreqSumFile = describe "readFreqSumFile" $ do+ (fsHeader, fsEntries) <- runIO . runSafeT $ do+ (fsHeader_, fsProd_) <- readFreqSumFile "testDat/example.freqsum"+ fsEntries_ <- purely P.fold list fsProd_+ return (fsHeader_, fsEntries_)+ it "should read the correct fs header" $+ fsHeader `shouldBe` mockDatFsHeader + it "should read the correct fs entries" $+ fsEntries `shouldBe` mockDatFsEntries ++testPrintFreqSumFile :: Spec+testPrintFreqSumFile = describe "printFreqSumFile" $ do+ let fn = "/tmp/freqSumWriteTest.txt"+ testDatFsProd = each mockDatFsEntries+ runIO . runSafeT . runEffect $ testDatFsProd >-> printFreqSumFile fn mockDatFsHeader+ (fsHeader, fsEntries) <- runIO . runSafeT $ do+ (fsHeader_, fsProd_) <- readFreqSumFile fn+ fsEntries_ <- purely P.fold list fsProd_+ return (fsHeader_, fsEntries_)+ it "should read the correct fs header after writing" $+ fsHeader `shouldBe` mockDatFsHeader + it "should read the correct fs entries after writing" $+ fsEntries `shouldBe` mockDatFsEntries ++mockDatFsHeader :: FreqSumHeader +mockDatFsHeader = FreqSumHeader names numbers+ where+ names = ["SAMPLE0", "SAMPLE1", "SAMPLE2", "SAMPLE3", "SAMPLE4"]+ numbers = [2, 2, 2, 1, 1]++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]]+
+ test/SequenceFormats/PileupSpec.hs view
@@ -0,0 +1,24 @@+module SequenceFormats.PileupSpec (spec) where++import SequenceFormats.Pileup (readPileupFromFile, PileupRow(..))+import SequenceFormats.Utils (Chrom(..))++import Control.Foldl (purely, list)+import qualified Pipes.Prelude as P+import Pipes.Safe (runSafeT)+import Test.Hspec++spec :: Spec+spec = testReadPileupFromFile++testReadPileupFromFile :: Spec+testReadPileupFromFile = describe "readPileupFromFile" $+ it "should read correct data from file" $ do+ let pProd = readPileupFromFile "testDat/example.pileup"+ runSafeT (purely P.fold list pProd) `shouldReturn` mockDatPentries++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"]]
+ test/SequenceFormats/RareAlleleHistogramSpec.hs view
@@ -0,0 +1,44 @@+{-# LANGUAGE OverloadedStrings #-}+module SequenceFormats.RareAlleleHistogramSpec (spec) where++import SequenceFormats.RareAlleleHistogram (RareAlleleHistogram(..), readHistogram, writeHistogramFile)++import qualified Data.Map as Map+import Test.Hspec++spec :: Spec+spec = do+ testReadHistogram+ testWriteHistogramFile++testReadHistogram :: Spec+testReadHistogram = describe "readHistogram" $+ it "should read correct data" $+ readHistogram "testDat/example.histogram.txt" `shouldReturn` mockHistogramDat++testWriteHistogramFile :: Spec+testWriteHistogramFile = describe "writeHistogramFile" $+ it "should read the correct histogram after writing" $ do+ let fn = "/tmp/histogramWriteTest.txt"+ writeHistogramFile fn mockHistogramDat+ readHistogram fn `shouldReturn`mockHistogramDat++mockHistogramDat :: RareAlleleHistogram+mockHistogramDat = RareAlleleHistogram names nVec 1 10 [] [] 1146826657 counts jnEstimates+ where+ names = ["EUR", "SEA", "SIB", "CHK", "SAM"]+ nVec = [66, 42, 44, 8, 28]+ counts = Map.fromList [+ ([1,0,0,0,0], 773148),+ ([0,1,0,0,0], 527207),+ ([0,0,1,0,0], 368640),+ ([0,0,0,0,1], 213918),+ ([2,0,0,0,0], 158795),+ ([0,2,0,0,0], 73241)]+ jnEstimates = Just $ Map.fromList [+ ([1,0,0,0,0], (6.7e-4, 3.6e-6)),+ ([0,1,0,0,0], (4.5e-4, 2.5e-6)),+ ([0,0,1,0,0], (3.2e-4, 2.4e-6)),+ ([0,0,0,0,1], (1.8e-4, 3.6e-6)),+ ([2,0,0,0,0], (1.3e-4, 1.4e-6)),+ ([0,2,0,0,0], (6.3e-5, 7.3e-7))]
+ test/SequenceFormats/UtilsSpec.hs view
@@ -0,0 +1,20 @@+{-# LANGUAGE OverloadedStrings #-}+module SequenceFormats.UtilsSpec (spec) where++import SequenceFormats.Utils (Chrom(..))++import Test.Hspec++spec :: Spec+spec = testChrom++testChrom :: Spec+testChrom = describe "Chrom" $ do+ specify "2 should be smaller than 10" $+ Chrom "2" < Chrom "10" `shouldBe` True+ specify "chr2 should be smaller than 10" $+ Chrom "chr2" < Chrom "10" `shouldBe` True+ specify "2 should be smaller than chr10" $+ Chrom "2" < Chrom "chr10" `shouldBe` True+ specify "chr2 should be smaller than chr10" $+ Chrom "chr2" < Chrom "chr10" `shouldBe` True
+ test/SequenceFormats/VCFSpec.hs view
@@ -0,0 +1,85 @@+{-# LANGUAGE OverloadedStrings #-}+module SequenceFormats.VCFSpec (spec) where++import Control.Foldl (list, purely)+import Pipes.Prelude (fold)+import Pipes.Safe (runSafeT)+import SequenceFormats.FreqSum (FreqSumEntry(..))+import SequenceFormats.Utils (Chrom(..))+import SequenceFormats.VCF (readVCFfromFile, getGenotypes, getDosages,+ isTransversionSnp, vcfToFreqSumEntry, isBiallelicSnp, VCFheader(..), VCFentry(..))+import Test.Hspec++spec :: Spec+spec = do+ testReadVCFfromFile+ testGetGenotypes+ testGetDosages+ testIsTransversionSnp+ testVcfToFreqsumEntry+ testIsBiallelicSnp++testReadVCFfromFile :: Spec+testReadVCFfromFile = describe "readVCFfromFile" $ do+ (vcfH, vcfRows) <- runIO . runSafeT $ do+ (vcfH_, vcfProd_) <- readVCFfromFile "testDat/example.vcf"+ vcfRows_ <- purely fold list vcfProd_+ return (vcfH_, vcfRows_)+ let vcfHc = vcfHeaderComments vcfH+ it "reads the correct header lines" $ do+ vcfHc !! 0 `shouldBe` "##fileformat=VCFv4.2"+ vcfHc !! 18 `shouldBe` "##bcftools_callCommand=call -c -v"+ it "reads the correct sample names" $+ vcfSampleNames vcfH `shouldBe` ["12880A", "12881A", "12883A", "12884A", "12885A"] + it "reads the correct vcf genotype rows" $ do+ vcfRows !! 0 `shouldBe` vcf1+ 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"]+ ["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"]+ ["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+testGetGenotypes = describe "getGenotypes" $ do+ it "should successfully read genotypes if GT format field is there" $+ getGenotypes vcf1 `shouldBe` Right ["0/0", "0/0", "0/1", "0/0", "0/0"]+ it "should yield Left err if GT format field isn't found" $+ getGenotypes (vcf1 {vcfFormatString=["PL"]}) `shouldBe` Left "GT format field not found"++testGetDosages :: Spec+testGetDosages = describe "getDosages" $ do+ it "should read correct dosages" $ do+ getDosages vcf1 `shouldBe` Right [Just 0, Just 0, Just 1, Just 0, Just 0]+ let vcf1' = vcf1 {vcfGenotypeInfo=[+ ["0/0", "0,3,37"], ["0/0", "0,6,67"], [".", "51,0,28"], ["0/0", "0,54,255"],+ ["0/0", "0,9,83"]]}+ getDosages vcf1' `shouldBe` Right [Just 0, Just 0, Nothing, Just 0, Just 0]++testIsTransversionSnp :: Spec+testIsTransversionSnp = describe "isTransversionSnp" $ do+ it "should reject triAllelic SNPs" $+ isTransversionSnp "A" ["C", "T"] `shouldBe` False+ it "should reject transitions" $ do+ isTransversionSnp "A" ["G"] `shouldBe` False+ isTransversionSnp "T" ["C"] `shouldBe` False+ it "should accept transversions" $ do+ isTransversionSnp "T" ["A"] `shouldBe` True+ isTransversionSnp "C" ["G"] `shouldBe` True++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])+ vcfToFreqSumEntry vcf1 `shouldBe` r++testIsBiallelicSnp :: Spec+testIsBiallelicSnp = describe "isBiallelicSnp" $ do+ it "should reject triAllelic" $+ isBiallelicSnp "A" ["C", "T"] `shouldBe` False+ it "should accept biallelic" $+ isBiallelicSnp "A" ["T"] `shouldBe` True
+ test/Spec.hs view
@@ -0,0 +1,1 @@+{-# OPTIONS_GHC -F -pgmF hspec-discover #-}