sequence-formats 1.1.4.1 → 1.3.1
raw patch · 15 files changed
+643/−216 lines, 15 filesdep +sequence-formatsdep +tastydep +tasty-hunitdep −MissingHdep −data-memocombinatorsdep −hsloggerdep ~attoparsecdep ~basedep ~bytestring
Dependencies added: sequence-formats, tasty, tasty-hunit
Dependencies removed: MissingH, data-memocombinators, hslogger, pipes-text, split, text, turtle
Dependency ranges changed: attoparsec, base, bytestring, containers, errors, exceptions, foldl, lens-family, pipes, pipes-attoparsec, pipes-bytestring, pipes-safe, transformers, vector
Files
- Changelog.md +22/−0
- sequence-formats.cabal +39/−21
- src-tests/SequenceFormats/Eigenstrat/Test.hs +79/−0
- src-tests/SequenceFormats/Fasta/Test.hs +17/−0
- src-tests/SequenceFormats/FreqSum/Test.hs +49/−0
- src-tests/SequenceFormats/RareAlleleHistogram/Test.hs +40/−0
- src-tests/SequenceFormats/Utils/Test.hs +12/−0
- src-tests/SequenceFormats/VCF/Test.hs +61/−0
- src-tests/tests.hs +48/−0
- src/SequenceFormats/Eigenstrat.hs +111/−63
- src/SequenceFormats/Fasta.hs +9/−10
- src/SequenceFormats/FreqSum.hs +23/−25
- src/SequenceFormats/RareAlleleHistogram.hs +42/−29
- src/SequenceFormats/Utils.hs +33/−18
- src/SequenceFormats/VCF.hs +58/−50
Changelog.md view
@@ -1,1 +1,23 @@ V 1.1.4.1: First entry in the Changelog. Added Haddock documentation to all modules and prepare for releasing on Hackage.++V 1.1.4.2: Exporting readVCFfromProd++V 1.1.5: Fixed VCF parser: Now breaks if lines end prematurely++V 1.1.6: VCF parser now allows for truncated VCF files with no Format and Genotypes (sites-only VCF files)++V 1.1.7: Added option to parse Bim file (slightly different layout to Eigenstrat Snp Format), and added genetic position and snpId to the EigenstratSnpEntry datastructure. This will cause breaking changes in code linking against this library.++V 1.1.8: Added more consumers for Bim and Eigenstrat Snp files++V 1.1.8.1: Added Eq class to EigenstratInd and Sex++V 1.1.8.2: Added Eq and Show classes to various FreqSum entities. Fixed writing function, added tests.++V 1.1.8.3: Added tests for Fasta import. Succeed.++V 1.2.0: Added tests for VCF, and several bugfixes. Now runs on LTS-14.1 with pipes-text as legacy dependency.++V 1.3.0: Removed pipes-text, text and turtle dependencies and some more. Restructured all datatypes to use Bytestring instead of text. ++V 1.3.1: Moved test suite outside of the main library into the test source directory. Cleaner setup.
sequence-formats.cabal view
@@ -1,6 +1,6 @@ cabal-version: >=1.10 name: sequence-formats-version: 1.1.4.1+version: 1.3.1 license: GPL-3 license-file: LICENSE maintainer: stephan.schiffels@mac.com@@ -26,23 +26,41 @@ default-language: Haskell2010 build-depends: base >=4.7 && <5,- containers >=0.5.10.2 && <0.6,- errors >=2.2.2 && <2.3,- text >=1.2.2.2 && <1.3,- attoparsec >=0.13.2.2 && <0.14,- pipes >=4.3.7 && <4.4,- pipes-text >=0.0.2.5 && <0.1,- transformers >=0.5.2.0 && <0.6,- turtle >=1.4.5 && <1.5,- bytestring >=0.10.8.2 && <0.11,- lens-family >=1.2.2 && <1.3,- pipes-bytestring >=2.1.6 && <2.2,- foldl >=1.3.7 && <1.4,- exceptions >=0.8.3 && <0.9,- pipes-safe >=2.2.6 && <2.3,- pipes-attoparsec >=0.5.1.5 && <0.6,- vector >=0.12.0.1 && <0.13,- data-memocombinators >=0.5.1 && <0.6,- split >=0.2.3.3 && <0.3,- MissingH >=1.4.0.1 && <1.5,- hslogger >=1.2.10 && <1.3+ containers >=0.6.0.1,+ errors >=2.3.0,+ attoparsec >=0.13.2.2,+ pipes >=4.3.11,+ transformers >=0.5.6.2,+ bytestring >=0.10.8.2,+ lens-family >=1.2.3,+ pipes-bytestring >=2.1.6,+ foldl >=1.4.5,+ exceptions >=0.10.2,+ pipes-safe >=2.3.1,+ pipes-attoparsec >=0.5.1.5,+ vector >=0.12.0.3++test-suite sequenceFormatTests+ type: exitcode-stdio-1.0+ main-is: tests.hs+ hs-source-dirs: src-tests+ other-modules:+ SequenceFormats.Eigenstrat.Test+ SequenceFormats.Fasta.Test+ SequenceFormats.FreqSum.Test+ SequenceFormats.RareAlleleHistogram.Test+ SequenceFormats.Utils.Test+ SequenceFormats.VCF.Test+ default-language: Haskell2010+ build-depends:+ base >=4.12.0.0,+ sequence-formats -any,+ foldl >=1.4.5,+ pipes >=4.3.11,+ pipes-safe >=2.3.1,+ tasty >=1.2.3,+ vector >=0.12.0.3,+ transformers >=0.5.6.2,+ tasty-hunit >=0.10.0.2,+ bytestring >=0.10.8.2,+ containers >=0.6.0.1
+ src-tests/SequenceFormats/Eigenstrat/Test.hs view
@@ -0,0 +1,79 @@+{-# 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 view
@@ -0,0 +1,17 @@+{-# 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 view
@@ -0,0 +1,49 @@+{-# 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 = "/Users/schiffels/repos/github/stschiff/sequence-formats/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 view
@@ -0,0 +1,40 @@+{-# 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 view
@@ -0,0 +1,12 @@+{-# 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 view
@@ -0,0 +1,61 @@+{-# 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 view
@@ -0,0 +1,48 @@+{-# 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
@@ -6,36 +6,45 @@ module SequenceFormats.Eigenstrat (EigenstratSnpEntry(..), EigenstratIndEntry(..), readEigenstratInd, GenoEntry(..), GenoLine, Sex(..), - readEigenstratSnpStdIn, readEigenstratSnpFile, readEigenstrat, writeEigenstrat) where+ readEigenstratSnpStdIn, readEigenstratSnpFile, readBimStdIn, readBimFile,+ readEigenstrat, writeEigenstrat, writeEigenstratIndFile, writeEigenstratSnp, writeBim, + writeEigenstratGeno) where -import SequenceFormats.Utils (consumeProducer, FormatException(..), Chrom(..))+import SequenceFormats.Utils (consumeProducer, SeqFormatException(..), Chrom(..), readFileProd) import Control.Applicative ((<|>))+import Control.Exception (throw) import Control.Monad (void, forM_)-import Control.Monad.Catch (MonadThrow, throwM)+import Control.Monad.Catch (MonadThrow) import Control.Monad.IO.Class (MonadIO, liftIO)-import qualified Data.Attoparsec.Text as A+import qualified Data.Attoparsec.ByteString.Char8 as A import Data.Char (isSpace) import Data.Vector (Vector, fromList, toList)-import qualified Data.Text as T-import qualified Data.Text.IO as T+import qualified Data.ByteString.Char8 as B import Pipes (Producer, Pipe, (>->), for, cat, yield, Consumer) import Pipes.Safe (MonadSafe)+import qualified Pipes.Safe.Prelude as PS import qualified Pipes.Prelude as P-import qualified Pipes.Text.IO as PT-import System.IO (withFile, IOMode(..))-import Turtle (format, w, d, (%), s)+import qualified Pipes.ByteString as PB+import System.IO (withFile, IOMode(..), Handle, hPutStrLn) -- |A datatype to represent a single genomic SNP. The constructor arguments are: -- Chromosome, Position, Reference Allele, Alternative Allele.-data EigenstratSnpEntry = EigenstratSnpEntry Chrom Int Char Char deriving (Eq, Show)+data EigenstratSnpEntry = EigenstratSnpEntry {+ snpChrom :: Chrom,+ snpPos :: Int,+ snpGeneticPos :: Double,+ snpId :: B.ByteString,+ snpRef :: Char,+ snpAlt :: Char+ } deriving (Eq, Show) -- |A datatype to represent a single individual. The constructor arguments are: -- Name, Sex and Population Name-data EigenstratIndEntry = EigenstratIndEntry T.Text Sex T.Text deriving (Show)+data EigenstratIndEntry = EigenstratIndEntry String Sex String deriving (Eq, Show) -- |A datatype to represent Sex in an Eigenstrat Individual file-data Sex = Male | Female | Unknown deriving (Show)+data Sex = Male | Female | Unknown deriving (Eq, Show) -- |A datatype to represent the genotype of an individual at a SNP. data GenoEntry = HomRef | Het | HomAlt | Missing deriving (Eq, Show)@@ -45,22 +54,27 @@ eigenstratSnpParser :: A.Parser EigenstratSnpEntry eigenstratSnpParser = do- A.skipMany A.space- void word- A.skipMany1 A.space- chrom <- word- A.skipMany1 A.space- void word- A.skipMany1 A.space- pos <- A.decimal- A.skipMany1 A.space- ref <- A.satisfy (A.inClass "ACTGN")- A.skipMany1 A.space- alt <- A.satisfy (A.inClass "ACTGX")+ snpId_ <- A.skipMany A.space >> word+ chrom <- A.skipMany1 A.space >> word+ geneticPos <- A.skipMany1 A.space >> A.double+ pos <- A.skipMany1 A.space >> A.decimal+ ref <- A.skipMany1 A.space >> A.satisfy (A.inClass "ACTGN")+ alt <- A.skipMany1 A.space >> A.satisfy (A.inClass "ACTGX") void A.endOfLine- return $ EigenstratSnpEntry (Chrom chrom) pos ref alt+ return $ EigenstratSnpEntry (Chrom (B.unpack chrom)) pos geneticPos snpId_ ref alt -word :: A.Parser T.Text+bimParser :: A.Parser EigenstratSnpEntry+bimParser = do+ chrom <- word+ snpId_ <- A.skipMany1 A.space >> word+ geneticPos <- A.skipMany1 A.space >> A.double+ pos <- A.skipMany1 A.space >> A.decimal+ ref <- A.skipMany1 A.space >> A.satisfy (A.inClass "ACTGN")+ alt <- A.skipMany1 A.space >> A.satisfy (A.inClass "ACTGX")+ 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@@ -72,7 +86,7 @@ A.skipMany1 A.space popName <- word void A.endOfLine- return $ EigenstratIndEntry name sex popName+ return $ EigenstratIndEntry (B.unpack name) sex (B.unpack popName) parseSex :: A.Parser Sex parseSex = parseMale <|> parseFemale <|> parseUnknown@@ -83,16 +97,16 @@ -- |Function to read an Eigenstrat individual file. Returns the Eigenstrat Individual Entries as list. readEigenstratInd :: (MonadIO m) => FilePath -> m [EigenstratIndEntry]-readEigenstratInd fn = do- liftIO . withFile fn ReadMode $ \handle -> do- P.toListM $ consumeProducer eigenstratIndParser (PT.fromHandle handle)+readEigenstratInd fn =+ liftIO . withFile fn ReadMode $ \handle ->+ P.toListM $ consumeProducer eigenstratIndParser (PB.fromHandle handle) eigenstratGenoParser :: A.Parser GenoLine eigenstratGenoParser = do line <- A.takeWhile1 isValidNum void A.endOfLine return . fromList $ do- l <- T.unpack line+ l <- B.unpack line case l of '0' -> return HomAlt '1' -> return Het@@ -104,12 +118,20 @@ -- |Function to read a Snp File from StdIn. Returns a Pipes-Producer over the EigenstratSnpEntries. readEigenstratSnpStdIn :: (MonadThrow m, MonadIO m) => Producer EigenstratSnpEntry m ()-readEigenstratSnpStdIn = consumeProducer eigenstratSnpParser PT.stdin+readEigenstratSnpStdIn = consumeProducer eigenstratSnpParser PB.stdin -- |Function to read a Snp File from a file. Returns a Pipes-Producer over the EigenstratSnpEntries. readEigenstratSnpFile :: (MonadSafe m) => FilePath -> Producer EigenstratSnpEntry m ()-readEigenstratSnpFile = consumeProducer eigenstratSnpParser . PT.readFile+readEigenstratSnpFile = consumeProducer eigenstratSnpParser . readFileProd +-- |Function to read a Bim File from StdIn. Returns a Pipes-Producer over the EigenstratSnpEntries.+readBimStdIn :: (MonadThrow m, MonadIO m) => Producer EigenstratSnpEntry m ()+readBimStdIn = consumeProducer bimParser PB.stdin++-- |Function to read a Bim File from a file. Returns a Pipes-Producer over the EigenstratSnpEntries.+readBimFile :: (MonadSafe m) => FilePath -> Producer EigenstratSnpEntry m ()+readBimFile = consumeProducer bimParser . readFileProd+ -- |Function to read a full Eigenstrat database from files. Returns a pair of the Eigenstrat Individual Entries, and a joint Producer over the snp entries and the genotypes. readEigenstrat :: (MonadSafe m) => FilePath -- ^The Genotype file -> FilePath -- ^The Snp File@@ -118,20 +140,68 @@ readEigenstrat genoFile snpFile indFile = do indEntries <- readEigenstratInd indFile let snpProd = readEigenstratSnpFile snpFile- genoProd = consumeProducer eigenstratGenoParser (PT.readFile genoFile) >-> + genoProd = consumeProducer eigenstratGenoParser (readFileProd genoFile) >-> validateEigenstratEntries (length indEntries) return (indEntries, P.zip snpProd genoProd) validateEigenstratEntries :: (MonadThrow m) => Int -> Pipe GenoLine GenoLine m ()-validateEigenstratEntries nr = for cat $ \line -> do+validateEigenstratEntries nr = for cat $ \line -> if length line /= nr then do- let msg = format ("inconsistent nr of genotypes ("%d%", but should be "%d%") in \- \genotype line "%w) (length line) nr line- throwM $ FormatException msg+ let msg = "inconsistent nr of genotypes (" <> show (length line) <> ", but should be " <> show nr <> ") in \+ \genotype line " <> show line+ throw (SeqFormatException msg) else yield line +-- |Function to write an Eigenstrat Ind file.+writeEigenstratIndFile :: (MonadIO m) => FilePath -> [EigenstratIndEntry] -> m ()+writeEigenstratIndFile f indEntries =+ liftIO . withFile f WriteMode $ \h ->+ forM_ indEntries $ \(EigenstratIndEntry name sex popName) ->+ hPutStrLn h $ name <> "\t" <> sexToStr sex <> "\t" <> popName+ where+ sexToStr sex = case sex of+ Male -> "M"+ Female -> "F"+ Unknown -> "U"++-- |Function to write an Eigenstrat Snp File. Returns a consumer expecting EigenstratSnpEntries.+writeEigenstratSnp :: (MonadIO m) => Handle -- ^The Eigenstrat Snp File Handle.+ -> Consumer EigenstratSnpEntry m () -- ^A consumer to read EigenstratSnpEntries+writeEigenstratSnp snpFileH =+ let snpOutTextConsumer = PB.toHandle snpFileH+ toTextPipe = P.map (\(EigenstratSnpEntry chrom pos gpos gid ref alt) ->+ let snpLine = B.intercalate "\t" [gid, B.pack (unChrom chrom), B.pack (show gpos),+ B.pack (show pos), B.singleton ref, B.singleton alt]+ in snpLine <> "\n")+ in toTextPipe >-> snpOutTextConsumer++-- |Function to write a Bim file. Returns a consumer expecting EigenstratSnpEntries.+writeBim :: (MonadIO m) => Handle -- ^The Eigenstrat Snp File handle.+ -> Consumer EigenstratSnpEntry m () -- ^A consumer to read EigenstratSnpEntries+writeBim snpFileH =+ let snpOutTextConsumer = PB.toHandle snpFileH+ toTextPipe = P.map (\(EigenstratSnpEntry chrom pos gpos gid ref alt) ->+ B.intercalate "\t" [B.pack (unChrom chrom), gid, B.pack (show gpos), B.pack (show pos), B.singleton ref, B.singleton alt])+ in toTextPipe >-> snpOutTextConsumer++-- |Function to write an Eigentrat Geno File. Returns a consumer expecting Eigenstrat Genolines.+writeEigenstratGeno :: (MonadIO m) => Handle -- ^The Genotype file handle+ -> Consumer GenoLine m () -- ^A consumer to read Genotype entries.+writeEigenstratGeno genoFileH =+ let genoOutTextConsumer = PB.toHandle genoFileH+ toTextPipe = P.map (\genoLine ->+ let genoLineStr = B.concat . map (B.pack . show . toEigenStratNum) . toList $ genoLine+ in genoLineStr <> "\n")+ in toTextPipe >-> genoOutTextConsumer+ where+ toEigenStratNum c = case c of+ HomRef -> 2 :: Int+ Het -> 1+ HomAlt -> 0+ Missing -> 9+ -- |Function to write an Eigenstrat Database. Returns a consumer expecting joint Snp- and Genotype lines. writeEigenstrat :: (MonadSafe m) => FilePath -- ^The Genotype file -> FilePath -- ^The Snp File@@ -140,29 +210,7 @@ -> Consumer (EigenstratSnpEntry, GenoLine) m () -- ^A consumer to read joint Snp/Genotype entries. writeEigenstrat genoFile snpFile indFile indEntries = do liftIO $ writeEigenstratIndFile indFile indEntries- let snpOutTextConsumer = PT.writeFile snpFile- genoOutTextConsumer = PT.writeFile genoFile- toTextPipe = P.map (\(EigenstratSnpEntry chrom pos ref alt, genoLine) ->- let n = format (s%"_"%d) (unChrom chrom) pos- snpLine = format (s%"\t"%s%"\t0\t"%d%"\t"%s%"\t"%s%"\n") n (unChrom chrom) pos- (T.singleton ref) (T.singleton alt)- genoLineStr = T.concat . map (format d . toEigenStratNum) . toList $ genoLine- in (snpLine, format (s%"\n") genoLineStr))- toTextPipe >-> P.tee (P.map fst >-> snpOutTextConsumer) >-> P.map snd >-> genoOutTextConsumer- where- toEigenStratNum c = case c of- HomRef -> 2 :: Int- Het -> 1- HomAlt -> 0- Missing -> 9+ let snpOutConsumer = PS.withFile snpFile WriteMode writeEigenstratSnp+ genoOutConsumer = PS.withFile genoFile WriteMode writeEigenstratGeno+ P.tee (P.map fst >-> snpOutConsumer) >-> P.map snd >-> genoOutConsumer -writeEigenstratIndFile :: (MonadIO m) => FilePath -> [EigenstratIndEntry] -> m ()-writeEigenstratIndFile f indEntries = do- liftIO . withFile f WriteMode $ \h -> do- forM_ indEntries $ \(EigenstratIndEntry name sex popName) -> do- liftIO . T.hPutStrLn h $ format (s%"\t"%s%"\t"%s) name (sexToStr sex) popName- where- sexToStr sex = case sex of- Male -> "M"- Female -> "F"- Unknown -> "U"
src/SequenceFormats/Fasta.hs view
@@ -14,16 +14,12 @@ import qualified Data.Attoparsec.ByteString.Char8 as A import qualified Data.ByteString.Char8 as B import Data.Char (isAlphaNum)-import qualified Data.Text as T import Lens.Family2 (view) import Pipes (Producer, next, (>->), runEffect) import Pipes.Attoparsec (parse) import qualified Pipes.ByteString as P import Pipes.Prelude (drain)-import System.IO (Handle)-import Turtle.Format (format, (%), s)-import Turtle.Prelude (err)-import Turtle.Line (unsafeTextToLine)+import System.IO (Handle, hPutStr, stderr) -- |A function to select out a specific chromosome from a Fasta File. Expects a file handle to the -- file and a chromosome. Note that by Chromosome I simply denote a fasta header line, as is the @@ -35,7 +31,7 @@ where go prod = do (chrom_, prod') <- readNextFastaEntry prod- err . unsafeTextToLine $ format ("found chromosome "%s) (unChrom chrom_)+ hPutStr stderr ("found chromosome " <> unChrom chrom_) if chrom_ == chrom then return (void prod') else do@@ -50,9 +46,12 @@ header <- case res of Nothing -> liftIO . throwIO $ AssertionFailed "Could not find chromosome. Fasta file exhausted." Just (Left e_) -> do- Right (chunk, _) <- next rest- let msg = show e_ ++ B.unpack chunk- liftIO . throwIO $ AssertionFailed ("Fasta header parsing error: " ++ msg)+ x <- next rest+ case x of+ (Right (chunk, _)) -> do+ let msg = show e_ ++ B.unpack chunk+ liftIO . throwIO $ AssertionFailed ("Fasta header parsing error: " ++ msg)+ _ -> error "should not happen" Just (Right h) -> return h return (header, view (P.break (==62)) rest >-> P.filter (\c -> c /= 10 && c /= 13)) -- '>' == 62, '\n' == 10, \r == 13@@ -64,4 +63,4 @@ A.skipSpace A.skipWhile (\c -> c /= '\n' && c /= '\r') A.endOfLine- return . Chrom . T.pack . B.unpack $ chrom+ return . Chrom . B.unpack $ chrom
src/SequenceFormats/FreqSum.hs view
@@ -13,31 +13,29 @@ import Control.Monad.Catch (MonadThrow, throwM) import Control.Monad.IO.Class (MonadIO, liftIO) import Control.Monad.Trans.State.Strict (runStateT)-import qualified Data.Attoparsec.Text as A+import qualified Data.Attoparsec.ByteString.Char8 as A import Data.Char (isAlphaNum, isSpace)-import Data.Text (Text, intercalate, singleton)-import Data.Text.IO (putStr, hPutStr)+import qualified Data.ByteString.Char8 as B import Pipes (Producer, (>->), Consumer) import Pipes.Attoparsec (parse, ParsingError(..)) import qualified Pipes.Prelude as P import Pipes.Safe (MonadSafe) import Pipes.Safe.Prelude (withFile)-import qualified Pipes.Text.IO as PT+import qualified Pipes.ByteString as PB import Prelude hiding (putStr) import System.IO (IOMode(..))-import Turtle (format, s, d, (%)) -- |A Datatype representing the Header data FreqSumHeader = FreqSumHeader {- fshNames :: [Text], -- ^A list of individual or group names+ fshNames :: [B.ByteString], -- ^A list of individual or group names fshCounts :: [Int] -- ^A list of haplotype counts per individual/group.-} deriving (Eq)+} deriving (Eq, Show) -freqSumHeaderToText :: FreqSumHeader -> Text+freqSumHeaderToText :: FreqSumHeader -> B.ByteString freqSumHeaderToText (FreqSumHeader names nCounts) =- format ("#CHROM\tPOS\tREF\tALT\t"%s%"\n") (intercalate "\t" tuples)+ "#CHROM\tPOS\tREF\tALT\t" <> B.intercalate "\t" tuples <> "\n" where- tuples = zipWith (\n c -> format (s%"("%d%")") n c) names nCounts+ tuples = zipWith (\n c -> n <> "(" <> B.pack (show c) <> ")") names nCounts -- |A Datatype to denote a single freqSum line data FreqSumEntry = FreqSumEntry {@@ -46,20 +44,19 @@ 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.-}+} deriving (Eq, Show) -- |This function converts a single freqSum entry to a printable freqSum line.-freqSumEntryToText :: FreqSumEntry -> Text+freqSumEntryToText :: FreqSumEntry -> B.ByteString freqSumEntryToText (FreqSumEntry chrom pos ref alt maybeCounts) =- format (s%"\t"%d%"\t"%s%"\t"%s%"\t"%s%"\n") (unChrom chrom) pos (singleton ref) (singleton alt) - countStr + B.intercalate "\t" [B.pack (unChrom chrom), B.pack (show pos), B.singleton ref, B.singleton alt, countStr] <> "\n" where- countStr = intercalate "\t" . map (format d . convertToNum) $ maybeCounts + countStr = B.intercalate "\t" . map (B.pack . show . convertToNum) $ maybeCounts convertToNum Nothing = -1 convertToNum (Just a) = a readFreqSumProd :: (MonadThrow m) =>- Producer Text m () -> m (FreqSumHeader, Producer FreqSumEntry m ())+ Producer B.ByteString m () -> m (FreqSumHeader, Producer FreqSumEntry m ()) readFreqSumProd prod = do (res, rest) <- runStateT (parse parseFreqSumHeader) prod header <- case res of@@ -70,11 +67,11 @@ -- |A function to read a freqsum file from StdIn. Returns a pair of a freqSum Header and a Producer over all lines. readFreqSumStdIn :: (MonadIO m, MonadThrow m) => m (FreqSumHeader, Producer FreqSumEntry m ())-readFreqSumStdIn = readFreqSumProd PT.stdin+readFreqSumStdIn = readFreqSumProd PB.stdin -- |A function to read a freqsum file from a file. Returns a pair of a freqSum Header and a Producer over all lines. readFreqSumFile :: (MonadSafe m) => FilePath -> m (FreqSumHeader, Producer FreqSumEntry m ())-readFreqSumFile = readFreqSumProd . PT.readFile+readFreqSumFile file = readFreqSumProd $ withFile file ReadMode PB.fromHandle parseFreqSumHeader :: A.Parser FreqSumHeader parseFreqSumHeader = do@@ -86,7 +83,7 @@ tuple = (,) <$> A.takeWhile (\c -> isAlphaNum c || c == '_' || c == '-') <* A.char '(' <*> A.decimal <* A.char ')' parseFreqSumEntry :: A.Parser FreqSumEntry-parseFreqSumEntry = FreqSumEntry <$> (Chrom <$> A.takeTill isSpace) <* A.skipSpace <*> A.decimal <*+parseFreqSumEntry = FreqSumEntry <$> (Chrom . B.unpack <$> A.takeTill isSpace) <* A.skipSpace <*> A.decimal <* A.skipSpace <*> base <* A.skipSpace <*> baseOrDot <* A.skipSpace <*> counts <* A.endOfLine where counts = (parseMissing <|> parseCount) `A.sepBy` A.char '\t'@@ -98,12 +95,13 @@ -- |A function to write freqSum data to StdOut. Expects the freqSum header as argument, and then returns a Consumer that accepts freqSum entries. printFreqSumStdOut :: (MonadIO m) => FreqSumHeader -> Consumer FreqSumEntry m () printFreqSumStdOut fsh = do- liftIO . putStr . freqSumHeaderToText $ fsh- P.map freqSumEntryToText >-> PT.stdout+ liftIO . B.putStr . freqSumHeaderToText $ fsh+ P.map freqSumEntryToText >-> PB.stdout -- |A function that writes a freqSum file. Expects the FilePath and the freqSum header as arguments, and then returns a Consumer that accepts freqSum entries. printFreqSumFile :: (MonadSafe m) => FilePath -> FreqSumHeader -> Consumer FreqSumEntry m ()-printFreqSumFile outFile fsh = do- outFileH <- withFile outFile WriteMode return- liftIO . hPutStr outFileH . freqSumHeaderToText $ fsh- P.map freqSumEntryToText >-> PT.toHandle outFileH+printFreqSumFile outFile fsh = withFile outFile WriteMode go+ where+ go h = do+ liftIO . B.hPutStr h . freqSumHeaderToText $ fsh+ P.map freqSumEntryToText >-> PB.toHandle h
src/SequenceFormats/RareAlleleHistogram.hs view
@@ -5,25 +5,28 @@ -} module SequenceFormats.RareAlleleHistogram (RareAlleleHistogram(..), readHistogramFromHandle,- SitePattern, readHistogram, showHistogram, showSitePattern) where+ SitePattern, readHistogram, writeHistogramStdOut, writeHistogramFile, showSitePattern) where +import SequenceFormats.Utils (SeqFormatException(..))+ import Control.Applicative (optional)-import Control.Error (Script, scriptIO, assertErr, throwE)+import Control.Error (assertErr)+import Control.Exception (throw)+import Control.Monad.IO.Class (MonadIO, liftIO) import Control.Monad.Trans.State.Strict (evalStateT)-import qualified Data.Attoparsec.Text as A+import qualified Data.Attoparsec.ByteString.Char8 as A import Data.Char (isAlphaNum) import Data.Int (Int64) import Data.List (intercalate, sortBy) import qualified Data.Map.Strict as Map-import qualified Data.Text as T+import qualified Data.ByteString.Char8 as B import Pipes.Attoparsec (parse)-import qualified Pipes.Text.IO as PT-import System.IO (Handle, openFile, IOMode(..), hClose)-import Turtle.Format ((%), w, format)+import qualified Pipes.ByteString as PB+import System.IO (Handle, IOMode(..), withFile) -- |A datatype to represent an Allele Sharing Histogram: data RareAlleleHistogram = RareAlleleHistogram {- raNames :: [T.Text], -- ^A list of branch names+ raNames :: [B.ByteString], -- ^A list of branch names raNVec :: [Int], -- ^A list of haploid sample sizes. raMinAf :: Int, -- ^The minimum allele count raMaxAf :: Int, -- ^The maximum allele count@@ -33,7 +36,7 @@ raTotalNrSites :: Int64, -- ^The total number of non-missing sites in the genome. raCounts :: Map.Map SitePattern Int64, -- ^The actual data, a dictionary from allele sharing patterns to observed numbers. raJackknifeEstimates :: Maybe (Map.Map SitePattern (Double, Double)) -- ^An optional dictionary that contains Jackknife estimates and standard deviations for each pattern frequency.-}+} deriving (Eq, Show) -- |A simple type synonym for the SitePattern, represented as a list of Integers that represents -- each pattern across the branches.@@ -41,49 +44,59 @@ -- |A simple function to convert a pattern into a String. showSitePattern :: SitePattern -> String-showSitePattern nVec = intercalate "," . map show $ nVec+showSitePattern = intercalate "," . map show -- |Function to convert a Rare Allele Histogram to text. Returns an error if attempting to print a -- histogram with non-standard settings. Many settings, such as minAf>1, are only meant for -- in-memory representations, but are not compatible with the file format itself.-showHistogram :: RareAlleleHistogram -> Either T.Text T.Text+showHistogram :: RareAlleleHistogram -> Either String B.ByteString showHistogram hist = do assertErr "can only print histogram with minAf=1 due to format-legacy" $ raMinAf hist == 1 assertErr "can only print histogram with no conditioning due to format-legacy" $ null (raConditionOn hist) assertErr "can only print histogram with no exclude pattern due to format-legacy" $ null (raExcludePatterns hist)- let head0 = T.concat ["NAMES=", T.intercalate "," . raNames $ hist]- head1 = T.concat ["N=", T.pack . intercalate "," . map show . raNVec $ hist]- head2 = T.concat ["MAX_M=", T.pack . show . raMaxAf $ hist]- head3 = T.concat ["TOTAL_SITES=", T.pack . show . raTotalNrSites $ hist]+ let head0 = "NAMES=" <> (B.intercalate "," . raNames $ hist)+ head1 = "N=" <> (B.pack . intercalate "," . map show . raNVec $ hist)+ head2 = "MAX_M=" <> (B.pack . show . raMaxAf $ hist)+ head3 = "TOTAL_SITES=" <> (B.pack . show . raTotalNrSites $ hist) body = do (k, v) <- sorted case raJackknifeEstimates hist of- Nothing -> [T.intercalate " " [T.pack . showSitePattern $ k, T.pack . show $ v]]+ Nothing -> [B.intercalate " " [B.pack . showSitePattern $ k, B.pack . show $ v]] Just jkHist -> do let Just (jkMean, jkSE) = k `Map.lookup` jkHist- return $ T.intercalate " " [T.pack . showSitePattern $ k, T.pack . show $ v,- T.pack . show $ jkMean, T.pack . show $ jkSE]- return $ T.unlines (head0:head1:head2:head3:body)+ return $ B.intercalate " " [B.pack . showSitePattern $ k, B.pack . show $ v,+ B.pack . show $ jkMean, B.pack . show $ jkSE]+ return $ B.unlines (head0:head1:head2:head3:body) where sorted = sortBy (\(_, v1) (_, v2) -> compare v2 v1) $ Map.toList (raCounts hist) +-- |Write a histogram to the stdout+writeHistogramStdOut :: (MonadIO m) => RareAlleleHistogram -> m ()+writeHistogramStdOut hist =+ case showHistogram hist of+ Left err -> throw (SeqFormatException err)+ Right outStr -> liftIO $ B.putStrLn outStr++-- |Write a histogram to a file+writeHistogramFile :: (MonadIO m) => FilePath -> RareAlleleHistogram -> m ()+writeHistogramFile outF hist =+ case showHistogram hist of+ Left err -> throw (SeqFormatException err)+ Right outStr -> liftIO $ B.writeFile outF outStr+ -- |Read a histogram from a FilePath-readHistogram :: FilePath -> Script RareAlleleHistogram-readHistogram path = do- h <- scriptIO $ openFile path ReadMode- hist <- readHistogramFromHandle h- scriptIO $ hClose h- return hist+readHistogram :: (MonadIO m) => FilePath -> m RareAlleleHistogram+readHistogram path = liftIO $ withFile path ReadMode readHistogramFromHandle -- |Read a histogram from a File Handle.-readHistogramFromHandle :: Handle -> Script RareAlleleHistogram+readHistogramFromHandle :: (MonadIO m) => Handle -> m RareAlleleHistogram readHistogramFromHandle handle = do- res <- evalStateT (parse parseHistogram) . PT.fromHandle $ handle+ res <- evalStateT (parse parseHistogram) . PB.fromHandle $ handle case res of- Nothing -> throwE "histogram file exhausted too early"- Just (Left err) -> throwE $ format ("Histogram parsing error: "%w) err+ Nothing -> throw (SeqFormatException "histogram file exhausted too early")+ Just (Left err) -> throw (SeqFormatException ("Histogram parsing error: " <> show err)) Just (Right hist) -> return hist parseHistogram :: A.Parser RareAlleleHistogram
src/SequenceFormats/Utils.hs view
@@ -3,20 +3,25 @@ -- |This module contains helper functions for file parsing. module SequenceFormats.Utils (liftParsingErrors,- consumeProducer,- FormatException(..),+ consumeProducer, readFileProd,+ SeqFormatException(..), Chrom(..)) where +import Control.Error (readErr) import Control.Exception (Exception) import Control.Monad.Catch (MonadThrow, throwM) import Control.Monad.Trans.Class (lift)-import qualified Data.Text as T-import qualified Data.Attoparsec.Text as A+import qualified Data.ByteString.Char8 as B+import qualified Data.Attoparsec.ByteString.Char8 as A import Pipes (Producer, next) import Pipes.Attoparsec (ParsingError(..), parsed)+import qualified Pipes.ByteString as PB+import qualified Pipes.Safe as PS+import qualified Pipes.Safe.Prelude as PS+import System.IO (IOMode(..)) -- |A wrapper datatype for Chromosome names.-newtype Chrom = Chrom {unChrom :: T.Text} deriving (Eq)+newtype Chrom = Chrom {unChrom :: String} deriving (Eq) -- |Show instance for Chrom instance Show Chrom where@@ -25,29 +30,39 @@ -- |Ord instance for Chrom instance Ord Chrom where compare (Chrom c1) (Chrom c2) = - let c1' = if T.take 3 c1 == "chr" then T.drop 3 c1 else c1- c2' = if T.take 3 c2 == "chr" then T.drop 3 c2 else c2- cn1 = read . T.unpack $ c1' :: Int- cn2 = read . T.unpack $ c2' :: Int- in cn1 `compare` cn2+ 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+ 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 FormatException = FormatException T.Text+data SeqFormatException = SeqFormatException String deriving Show -instance Exception FormatException+instance Exception SeqFormatException -- |A function to help with reporting parsing errors to stderr. Returns a clean Producer over the -- parsed datatype. liftParsingErrors :: (MonadThrow m) =>- Either (ParsingError, Producer T.Text m r) () -> Producer a m ()+ Either (ParsingError, Producer B.ByteString m r) () -> Producer a m () liftParsingErrors res = case res of- Left (ParsingError cont msg, restProd) -> do- Right (chunk, _) <- lift $ next restProd- let msg' = msg ++ " Error occurred while trying to parse this chunk: " ++ show chunk- throwM (ParsingError cont msg')+ Left (ParsingError _ msg, restProd) -> do+ x <- lift $ next restProd+ case x of+ Right (chunk, _) -> do+ let msg' = "Error while parsing: " <> msg <> ". Error occurred when trying to parse this chunk: " ++ show chunk+ throwM $ SeqFormatException msg'+ Left _ -> error "should not happen" Right () -> return () -- |A helper function to parse a text producer, properly reporting all errors to stderr.-consumeProducer :: (MonadThrow m) => A.Parser a -> Producer T.Text m () -> Producer a m ()+consumeProducer :: (MonadThrow m) => A.Parser a -> Producer B.ByteString m () -> Producer a m () consumeProducer parser prod = parsed parser prod >>= liftParsingErrors++readFileProd :: (PS.MonadSafe m) => FilePath -> Producer B.ByteString m ()+readFileProd f = PS.withFile f ReadMode PB.fromHandle+
src/SequenceFormats/VCF.hs view
@@ -8,33 +8,34 @@ VCFentry(..), readVCFfromStdIn, readVCFfromFile,+ readVCFfromProd, getGenotypes, getDosages, isTransversionSnp, vcfToFreqSumEntry, isBiallelicSnp) where -import SequenceFormats.Utils (consumeProducer, Chrom(..))+import SequenceFormats.Utils (consumeProducer, Chrom(..), readFileProd, SeqFormatException(..)) import SequenceFormats.FreqSum (FreqSumEntry(..)) -import Control.Applicative ((<|>), empty)+import Control.Applicative ((<|>)) import Control.Error (headErr, assertErr) import Control.Monad (void) import Control.Monad.Catch (MonadThrow, throwM) import Control.Monad.Trans.State.Strict (runStateT) import Control.Monad.IO.Class (MonadIO)-import qualified Data.Attoparsec.Text as A+import qualified Data.Attoparsec.ByteString.Char8 as A import Data.Char (isSpace)-import qualified Data.Text as T+import qualified Data.ByteString.Char8 as B import Pipes (Producer)-import Pipes.Attoparsec (parse, ParsingError(..))+import Pipes.Attoparsec (parse) import Pipes.Safe (MonadSafe)-import qualified Pipes.Text.IO as PT+import qualified Pipes.ByteString as PB -- |A datatype to represent the VCF Header. Most comments are simply parsed as entire lines, but the very last comment line, containing the sample names, is separated out data VCFheader = VCFheader {- vcfHeaderComments :: [T.Text], -- ^A list of containing all comments starting with a single '#'- vcfSampleNames :: [T.Text] -- ^The list of sample names parsed from the last comment line + vcfHeaderComments :: [String], -- ^A list of containing all comments starting with a single '#'+ vcfSampleNames :: [String] -- ^The list of sample names parsed from the last comment line -- starting with '##' } deriving (Show) @@ -42,70 +43,77 @@ data VCFentry = VCFentry { vcfChrom :: Chrom, -- ^The chromosome vcfPos :: Int, -- ^The position- vcfId :: Maybe T.Text, -- ^The SNP ID if non-missing- vcfRef :: T.Text, -- ^ The reference allele (supports also multi-character alleles for Indels)- vcfAlt :: [T.Text], -- ^The alternative alleles, each one possible of multiple characters + vcfId :: Maybe B.ByteString, -- ^The SNP ID if non-missing+ vcfRef :: B.ByteString, -- ^ The reference allele (supports also multi-character alleles for Indels)+ vcfAlt :: [B.ByteString], -- ^The alternative alleles, each one possible of multiple characters vcfQual :: Double, -- ^The quality value- vcfFilter :: Maybe T.Text, -- ^The Filter value, if non-missing.- vcfInfo :: [T.Text], -- ^A list of Info fields- vcfFormatString :: [T.Text], -- ^A list of format tags- vcfGenotypeInfo :: [[T.Text]] -- ^A list of format fields for each sapmle.-} deriving (Show)+ vcfFilter :: Maybe B.ByteString, -- ^The Filter value, if non-missing.+ vcfInfo :: [B.ByteString], -- ^A list of Info fields+ vcfFormatString :: [B.ByteString], -- ^A list of format tags+ vcfGenotypeInfo :: [[B.ByteString]] -- ^A list of format fields for each sample.+} deriving (Show, Eq) +-- |reads a VCFheader and VCFentries from a text producer. readVCFfromProd :: (MonadThrow m) =>- Producer T.Text m () -> m (VCFheader, Producer VCFentry m ())+ Producer B.ByteString m () -> m (VCFheader, Producer VCFentry m ()) readVCFfromProd prod = do (res, rest) <- runStateT (parse vcfHeaderParser) prod header <- case res of- Nothing -> throwM $ ParsingError [] "freqSum file exhausted"- Just (Left e) -> throwM e+ Nothing -> throwM $ SeqFormatException "freqSum file exhausted"+ Just (Left e) -> throwM (SeqFormatException (show e)) Just (Right h) -> return h return (header, consumeProducer vcfEntryParser rest) -- |Reading a VCF from StdIn. Returns a VCFHeader and a Producer over VCFentries. readVCFfromStdIn :: (MonadIO m, MonadThrow m) => m (VCFheader, Producer VCFentry m ())-readVCFfromStdIn = readVCFfromProd PT.stdin+readVCFfromStdIn = readVCFfromProd PB.stdin -- |Reading a VCF from a file. Returns a VCFHeader and a Producer over VCFentries. readVCFfromFile :: (MonadSafe m) => FilePath -> m (VCFheader, Producer VCFentry m ())-readVCFfromFile = readVCFfromProd . PT.readFile+readVCFfromFile = readVCFfromProd . readFileProd vcfHeaderParser :: A.Parser VCFheader vcfHeaderParser = VCFheader <$> A.many1' doubleCommentLine <*> singleCommentLine where doubleCommentLine = do c1 <- A.string "##"- s_ <- A.takeWhile1 (not . A.isEndOfLine)+ s_ <- A.takeWhile1 (/='\n') A.endOfLine- return $ T.append c1 s_+ return . B.unpack $ c1 <> s_ singleCommentLine = do void $ A.char '#'- s_ <- A.takeWhile1 (not . A.isEndOfLine)+ s_ <- A.takeWhile1 (/='\n') A.endOfLine- let fields = T.splitOn "\t" s_- return . drop 9 $ fields+ let fields = B.splitWith (=='\t') s_+ return . drop 9 . map B.unpack $ fields vcfEntryParser :: A.Parser VCFentry-vcfEntryParser = VCFentry <$> (Chrom <$> word) <* sp <*> A.decimal <* sp <*> parseId <* sp <*>- word <* sp <*> parseAlternativeAlleles <* sp <*> A.double <* sp <*> parseFilter <* sp <*> - parseInfoFields <* sp <*> parseFormatStrings <* sp <*> parseGenotypeInfos <* A.endOfLine+vcfEntryParser = vcfEntryParserFull <|> vcfEntryParserTruncated where- word = A.takeTill A.isHorizontalSpace- sp = A.satisfy A.isHorizontalSpace- parseId = parseDot <|> (Just <$> word)- parseDot = A.char '.' *> empty- parseAlternativeAlleles = parseDot <|> (parseAllele `A.sepBy1` A.char ',')- parseAllele = A.takeTill (\c -> c == ',' || A.isHorizontalSpace c)- parseFilter = parseDot <|> (Just <$> word)- parseInfoFields = parseDot <|> (parseInfoField `A.sepBy1` A.char ';')- parseInfoField = A.takeTill (\c -> c == ';' || A.isHorizontalSpace c)- parseFormatStrings = A.takeTill (\c -> c == ':' || A.isHorizontalSpace c) `A.sepBy1` A.char ':'- parseGenotypeInfos = parseGenotype `A.sepBy1` (A.satisfy A.isHorizontalSpace)- parseGenotype = parseGenoField `A.sepBy1` (A.char ':')+ vcfEntryParserFull = VCFentry <$> (Chrom . B.unpack <$> word) <* sp <*> A.decimal <* sp <*> parseId <*+ sp <*> word <* sp <*> parseAlternativeAlleles <* sp <*> A.double <* sp <*> parseFilter <* + sp <*> parseInfoFields <* sp <*> parseFormatStrings <* sp <*> parseGenotypeInfos <* + A.endOfLine+ 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 '.'+ parseAlternativeAlleles = (parseDot *> pure []) <|> (parseAllele `A.sepBy1` A.char ',')+ parseAllele = A.takeTill (\c -> c == ',' || isSpace c)+ parseFilter = (parseDot *> pure Nothing) <|> (Just <$> word)+ parseInfoFields = (parseDot *> pure []) <|> (parseInfoField `A.sepBy1` A.char ';')+ parseInfoField = A.takeTill (\c -> c == ';' || isSpace c)+ parseFormatStrings = parseFormatString `A.sepBy1` A.char ':'+ parseFormatString = A.takeTill (\c -> c == ':' || isSpace c)+ parseGenotypeInfos = parseGenotype `A.sepBy1` sp+ parseGenotype = parseGenoField `A.sepBy1` A.char ':' parseGenoField = A.takeTill (\c -> c == ':' || isSpace c) -- |returns True if the SNP is biallelic.-isBiallelicSnp :: T.Text -> [T.Text] -> Bool+isBiallelicSnp :: B.ByteString -> [B.ByteString] -> Bool isBiallelicSnp ref alt = validRef && validAlt where validRef = (ref `elem` ["A", "C", "G", "T"])@@ -114,7 +122,7 @@ _ -> False -- |returns True if the SNp is a biallelic Transversion SNP (i.e. one of G/T, G/C, A/T, A/C)-isTransversionSnp :: T.Text -> [T.Text] -> Bool+isTransversionSnp :: B.ByteString -> [B.ByteString] -> Bool isTransversionSnp ref alt = case alt of [alt'] -> isBiallelicSnp ref alt && (not $ isTransition ref alt')@@ -124,7 +132,7 @@ ((r == "C") && (a == "T")) || ((r == "T") && (a == "C")) -- |Extracts the genotype fields (for each sapmle) from a VCF entry-getGenotypes :: VCFentry -> Either String [T.Text]+getGenotypes :: VCFentry -> Either String [B.ByteString] getGenotypes vcfEntry = do gtIndex <- fmap fst . headErr "GT format field not found" . filter ((=="GT") . snd) . zip [0..] . vcfFormatString $ vcfEntry@@ -136,21 +144,21 @@ genotypes <- getGenotypes vcfEntry let dosages = do gen <- genotypes- if '.' `elem` (T.unpack gen) then+ if '.' `elem` (B.unpack gen) then return Nothing else- return . Just $ T.count "1" gen+ return . Just $ B.count '1' gen return dosages -- |Converts a VCFentry to the simpler FreqSum format (returns a Left Error if it fails.) vcfToFreqSumEntry :: VCFentry -> Either String FreqSumEntry vcfToFreqSumEntry vcfEntry = do dosages <- getDosages vcfEntry- assertErr "multi-site reference allele" $ T.length (vcfRef vcfEntry) == 1+ assertErr "multi-site reference allele" $ B.length (vcfRef vcfEntry) == 1 assertErr "need exactly one alternative allele" $ length (vcfAlt vcfEntry) == 1- assertErr "multi-site alternative allele" $ T.length (head . vcfAlt $ vcfEntry) == 1- let ref = T.head (vcfRef vcfEntry)- let alt = T.head . head . vcfAlt $ vcfEntry+ assertErr "multi-site alternative allele" $ B.length (head . vcfAlt $ vcfEntry) == 1+ let ref = B.head (vcfRef vcfEntry)+ 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