hyraxAbif (empty) → 0.2.3.2
raw patch · 23 files changed
+1854/−0 lines, 23 filesdep +basedep +binarydep +bytestringsetup-changedbinary-added
Dependencies added: base, binary, bytestring, directory, filepath, hedgehog, hscolour, hyraxAbif, pretty-show, protolude, text
Files
- LICENSE +30/−0
- README.md +265/−0
- Setup.hs +2/−0
- app/Main.hs +74/−0
- docs/eg_acag_acgt.png binary
- docs/eg_acag_acgt03.png binary
- docs/eg_acag_acgt_mix.png binary
- docs/eg_actg.png binary
- docs/eg_multi_mix.png binary
- hyraxAbif.cabal +85/−0
- src/Examples.hs +10/−0
- src/Examples/AddComment.hs +29/−0
- src/Examples/ReadAb1.hs +28/−0
- src/Examples/RemoveComments.hs +34/−0
- src/Hyrax/Abif.hs +124/−0
- src/Hyrax/Abif/Fasta.hs +50/−0
- src/Hyrax/Abif/Generate.hs +320/−0
- src/Hyrax/Abif/Read.hs +249/−0
- src/Hyrax/Abif/Write.hs +334/−0
- test/AbifTests.hs +119/−0
- test/FastaTests.hs +27/−0
- test/Generators.hs +53/−0
- test/Tests.hs +21/−0
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright HyraxBio (c) 2018++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++ * Redistributions in binary form must reproduce the above+ copyright notice, this list of conditions and the following+ disclaimer in the documentation and/or other materials provided+ with the distribution.++ * Neither the name of Author name here nor the names of other+ contributors may be used to endorse or promote products derived+ from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ README.md view
@@ -0,0 +1,265 @@+# HyraxBio AB1 parser, writer and generator (beta 0.2)++This project contains++ - Modules for parsing, generating or manipulating AB1 files.+ - Support for generating a minimal AB1 file from a FASTA input file+ - A simple terminal app to perform these operations++See++ - https://hackage.haskell.org/package/hyraxAbif for the hackage documentation+ - http://www6.appliedbiosystems.com/support/software_community/ABIF_File_Format.pdf for a high level overview of the AB1 file format.++# Building++Build with one of++ - `stack build` or (`make build`)+ - `cabal new-build`+++# Terminal app++Run with++ - `stack exec hyraxAbif-exe -- -- dump` if you are using stack+ - `cabal new-run hyraxAbif-exe dump` if you are using cabal 2.x++## Dump AB1++To dump an existing AB1 run++ `hyraxAbif-exe dump example.ab1`++This will output the structure of the AB1 like this++```+Header { hName = "ABIF" , hVersion = 101 }+Directory+ { dTagName = "tdir"+ , dTagNum = 1+ , dElemTypeCode = 1023+ , dElemTypeDesc = "root"+ , dElemType = ElemRoot+ , dElemSize = 28+ , dElemNum = 13+ , dDataSize = 364+ , dDataOffset = 61980+ , dData = ""+ , dDataDebug = []+ }+[ Directory+ { dTagName = "DATA"+ , dTagNum = 9+ , dElemTypeCode = 4+ , dElemTypeDesc = "short"+ , dElemType = ElemShort+ , dElemSize = 2+ , dElemNum = 7440+ , dDataSize = 14880+ , dDataOffset = 128+ , dData = ""+ , dDataDebug = []+ }+ +.+.+.++DATA {short} tagNum=9 size=2 count=7440 offset=128 []+DATA {short} tagNum=10 size=2 count=7440 offset=15008 []+DATA {short} tagNum=11 size=2 count=7440 offset=29888 []+DATA {short} tagNum=12 size=2 count=7440 offset=44768 []+FWO_ {char} tagNum=1 size=1 count=4 offset=1195463747 ["GATC"]+LANE {short} tagNum=1 size=2 count=1 offset=65536 ["1"]+PBAS {char} tagNum=1 size=1 count=744 offset=59648 ["GGGGGCAACTAAAGGAAGCTCTATTAGATACAGGAGCAGATGATACAGTATTAGAAGAAATGAGTTTGCCAGGAAGATGGAAACCAAAAATGATAGGGGGAATTGGAGGTTTTATCAAAGTAAGACAGTATGATCAGATACTCATAGAAATCTGTGGACATAAAGCTATAGGTACAGTATTAGTAGGACCTACACCTGTCAACATAATTGGAAGAAATCTGTTGACTCAGATTGGTTGCACTTTAAATTTTCCCATTAGCCCTATTGAGACTGTACCAGTAAAATTAAAGCCAGGAATGGATGGCCCAAAAGTTAAACAATGGCCATTGACAGAAGAAAAAATAAAAGCATTAGTAGAAATTTGTACAGAGATGGAAAAGGAAGGGAAAATTTCAAAAATTGGGCCTGAAAATCCATACAATACTCCAGTATTTGCCATAAAGAAAAAAGACAGTACTAAATGGAGAAAATTAGTAGATTTCAGAGAACTTAATAAGAGAACTCAAGACTTCTGGGAAGTTCAATTAGGAATACCACATCCCGCAGGGTTAAAAAAGAAAAAATCAGTAACAGTACTGGATGTGGGTGATGCATATTTTTCAGTTCCCTTAGATGAAGACTTCAGGAAGTATACTGCATTTACCATACCTAGTATAAACAATGAGACACCAGGGATTAGATATCAGTACAATGTGCTTCCACAGGGATGGAAAGGATCACCAGCAATATTCCAAAGTAGCATGA"]+PDMF {pString} tagNum=1 size=1 count=23 offset=60392 ["KB_3500_POP7_BDTv3.mob"]+PDMF {pString} tagNum=2 size=1 count=23 offset=60415 ["KB_3500_POP7_BDTv3.mob"]+PLOC {short} tagNum=1 size=2 count=744 offset=60438 []+S/N% {short} tagNum=1 size=2 count=4 offset=61926 []+SMPL {pString} tagNum=1 size=1 count=10 offset=61934 ["S17-SeqF1"]+CMNT {pString} tagNum=1 size=1 count=1 offset=61944 ["Generated by HyraxBio AB1 generator"]+```++The data is output twice. The first section is the detail, the second is the summary.++Selected data types have the "debug data" element populated. e.g. the PBAS (FASTA)++## Generate minimal AB1s from FASTAs++To create an AB1 run++ `hyraxAbif-exe gen "./pathContainingFastas" "./pathForOutputAb1s"`++This will create an AB1 per input FASTA++### Input FASTA format++Each input data should have the following format++```+> weight+read+> weight+read+```++ - The **weight** is a numeric value between 0 and 1 that specifies the weight of the current read. No other header/name is allowed++ - The **read** is the set of input nucleotides, IUPAC ambiguity codes are supported (MRWSYKVHDBNX). A read can be single or multi-line+++### Weighted reads++ - The weigh of a read specifies the intensity of the peak from 0 to 1. + - Weights for each position are added to a maximum of 1 per nucleotide+ - You can use `_` as a "blank" nucleotide, in which only the nucleotides from other reads will be considered+++For example+++ ```+ > 0.5+ ACG+ > 0.3+ AAAA+ > 1+ __AC+ ```+Results in the following weighted nucleotide per position++ * 0: `A` (0.5 + 0.3)+ * 1: `C` (0.5), `A` (0.3)+ * 2: `G` (0.5), `A` (0.3 + 1 = 1)+ * 3: `A` (0.3), `C` (1)+++*Note that the reads do not need to be the same length.*++---++#### Example FASTA - single file++***eg1.fasta***+```+> 1+ACTG+```+++++Here there is a single FASTA with a single read with a weigh of 1 (100%). The chromatogram for this AB1 shows perfect traces for the input `ACTG` nucleotides++---++#### Example FASTA - two FASTA files++***eg1.fasta***+```+> 1+ACAG+```++***eg2.fasta***+```+> 1+ACTG+```++++Two input FASTA files both with a weigh of 1. You can see in the second trace that the third nucleotide is a `T` (the trace is green). Exactly what the base-calling software (phred & recall etc) decide to call the base as depends on your settings and software choices.++---++#### Example FASTA - two FASTA files with different weights++***eg1.fasta***+```+> 1+ACAG+```++***eg2.fasta***+```+> 0.3+ACTG+```++++Here the second fasta has a weight of 0.3 and you can see the traces are 30% of the height of the top ones.+++---++#### Example FASTA - single FASTA with a mix++***eg1.fasta***+```+> 1+ACAG+> 0.3+ACTG+```++++The single input FASTA has an `AT` mix at the third nucleotide. The first read has a weight of 1 and the second a weight of 0.3.+Notice that the maximum weight is 1, e.g. the first `A` has the same intensity as the second even though the first one has the reads weighted both 1 and 0.3++---+++#### Example FASTA - Multiple mixes++***eg1.fasta***+```+> 1+ACAG+> 0.3+_GT+> 0.2+_G+```++++---++# Using the modules++ - Hyrax.Abif: The core AB1 types+ - Hyrax.Abif.Fasta: A simple FASTA parser used when generating AB1s+ - Hyrax.Abif.Read: Module for parsing an existing AB1+ - Hyrax.Abif.Write: Module for writing a new AB1 file+ - Hyrax.Abif.Generate: Module for generating a minimal AB1 from a given FASTA input++For a detailed overview of the code see *TODO* and the haddock documentation *TODO*++For now the terminal app (Main.hs) serves as an example and the best starting point to understand the code+++## E.g. Add a comment to an existing AB1 file+++```+import qualified Hyrax.Abif as H+import qualified Hyrax.Abif.Read as H+import qualified Hyrax.Abif.Write as H++addComment :: IO ()+addComment = do+ abif' <- H.readAbif "example.ab1"++ case abif' of+ Left e -> putText $ "error reading ABIF: " <> e+ Right abif -> do+ let modified = H.addDirectory abif $ H.mkComment "new comment"+ H.writeAbif "example.modified.ab1" modified+```++For additional examples see the Examples directory
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ app/Main.hs view
@@ -0,0 +1,74 @@+{-# LANGUAGE NoImplicitPrelude #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE LambdaCase #-}++module Main where++import Protolude+import qualified Data.Text as Txt+import qualified Data.ByteString.Lazy as BSL+import Text.Show.Pretty (ppShow)+import qualified System.Environment as Env+import qualified Language.Haskell.HsColour as Clr+import qualified Language.Haskell.HsColour.Colourise as Clr++import qualified Hyrax.Abif as H+import qualified Hyrax.Abif.Read as H+import qualified Hyrax.Abif.Generate as H++main :: IO ()+main = do+ args <- Env.getArgs ++ case headMay args of+ Just "dump" -> runDump+ Just "gen" -> runGenerateAb1+ _ -> putText "unknown args, expecting dump/gen"+ ++runGenerateAb1 :: IO ()+runGenerateAb1 =+ Env.getArgs >>= \case+ [_, source, dest] -> H.generateAb1s source dest+ _ -> putText "Expecting `source dest` args"+ ++runDump :: IO ()+runDump =+ Env.getArgs >>= \case+ [_, path] -> do+ file <- BSL.readFile path+ case H.getAbif file of+ Left e -> putText e+ Right (H.Abif hdr root dirs) -> do+ let debugged = H.clear . H.getDebug <$> dirs + --colourPrint debugged + colourPrint hdr+ colourPrint . H.clear $ root+ colourPrint debugged++ putText . Txt.intercalate "\n" $+ (\d -> H.dTagName d+ <> " {" <> H.dElemTypeDesc d <> "} tagNum="+ <> show (H.dTagNum d)+ <> " size=" <> show (H.dElemSize d)+ <> " count=" <> show (H.dElemNum d)+ <> " offset=" <> show (H.dDataOffset d)+ <> " "+ <> show (H.dDataDebug d)+ ) <$> debugged+ _ ->+ putText "Expecting path to ab1"+++myColourPrefs :: Clr.ColourPrefs+myColourPrefs = Clr.defaultColourPrefs { Clr.conop = [Clr.Foreground Clr.Yellow]+ , Clr.conid = [Clr.Foreground Clr.Yellow, Clr.Bold]+ , Clr.string = [Clr.Foreground $ Clr.Rgb 29 193 57]+ , Clr.char = [Clr.Foreground Clr.Cyan]+ , Clr.number = [Clr.Foreground $ Clr.Rgb 202 170 236]+ , Clr.keyglyph = [Clr.Foreground Clr.Yellow]+ }+ +colourPrint :: (Show a) => a -> IO ()+colourPrint = putStrLn . Clr.hscolour Clr.TTY myColourPrefs False False [] False . ppShow
+ docs/eg_acag_acgt.png view
binary file changed (absent → 5964 bytes)
+ docs/eg_acag_acgt03.png view
binary file changed (absent → 5789 bytes)
+ docs/eg_acag_acgt_mix.png view
binary file changed (absent → 4817 bytes)
+ docs/eg_actg.png view
binary file changed (absent → 3576 bytes)
+ docs/eg_multi_mix.png view
binary file changed (absent → 5050 bytes)
+ hyraxAbif.cabal view
@@ -0,0 +1,85 @@+cabal-version: 2.2+name: hyraxAbif+version: 0.2.3.2+synopsis: Modules for parsing, generating and manipulating AB1 files.+homepage: https://github.com/hyraxbio/hyraxAbif/#readme+license: BSD-3-Clause+license-file: LICENSE+author: HyraxBio+maintainer: andre@hyraxbio.co.za, andre@andrevdm.com+copyright: 2018 HyraxBio+category: Bioinformatics+build-type: Simple+extra-source-files: README.md+extra-doc-files: docs/*.png+description: This library provides functionality for parsing, modifying, writing and generating ABIF files+ <<docs/eg_multi_mix.png>>+ .+ Any AB1 file confirming to the standard at <http://www6.appliedbiosystems.com/support/software_community/ABIF_File_Format.pdf> should be supported.+ . + This library also support generating a minimal ABIF file from a FASTA input sequence.+ .+ A basic terminal application that can dump and generate AB1s is included. See <https://github.com/hyraxbio/hyraxAbif/blob/master/app/Main.hs>+ .+ See+ .+ * "Hyrax.Abif.Generate" for generate ABIF files from FASTA inputs+ * "Hyrax.Abif.Read" for parsing ABIF files+ * "Hyrax.Abif.Write" for creating/updating ABIF files+ * "Hyrax.Abif" for the core types+ * "Examples" for examples++library+ ghc-options: -Wall -Wincomplete-uni-patterns -Wincomplete-record-updates -Wimplicit-prelude -Wcompat -Wredundant-constraints -Wnoncanonical-monad-instances+ hs-source-dirs: src+ exposed-modules: Hyrax.Abif+ , Hyrax.Abif.Read+ , Hyrax.Abif.Write+ , Hyrax.Abif.Generate+ , Hyrax.Abif.Fasta+ , Examples+ , Examples.AddComment+ , Examples.RemoveComments+ , Examples.ReadAb1+ build-depends: base >= 4.9.1.0 && < 5+ , protolude >= 0.2.2 && < 0.2.3+ , text >= 1.2.3.0 && < 1.2.4.0+ , bytestring >= 0.10.8.2 && < 0.10.9.0+ , binary >= 0.8.5.1 && < 0.8.6.0+ , directory >= 1.3.0.2 && < 1.3.2.0+ , filepath >= 1.4.1.2 && < 1.4.2.0+ default-language: Haskell2010++executable hyraxAbif-exe+ hs-source-dirs: app+ main-is: Main.hs+ ghc-options: -threaded -rtsopts -with-rtsopts=-N -Wall -Wincomplete-uni-patterns -Wincomplete-record-updates -Wimplicit-prelude -Wcompat -Wredundant-constraints -Wnoncanonical-monad-instances+ build-depends: base >= 4.9.1.0 && < 5+ , hyraxAbif+ , protolude >= 0.2.2 && < 0.2.3+ , text >= 1.2.3.0 && < 1.2.4.0+ , bytestring >= 0.10.8.2 && < 0.10.9.0+ , pretty-show >= 1.6.16 && < 1.7.0+ , hscolour >= 1.24.4 && < 1.25.0+ default-language: Haskell2010++test-suite hyraxAbif-test+ type: exitcode-stdio-1.0+ hs-source-dirs: test+ main-is: Tests.hs+ build-depends: base >= 4.9.1.0 && < 5+ , hyraxAbif+ , protolude >= 0.2.2 && < 0.2.3+ , text >= 1.2.3.0 && < 1.2.4.0+ , bytestring >= 0.10.8.2 && < 0.10.9.0+ , binary >= 0.8.5.1 && < 0.8.6.0+ , hedgehog >= 0.5.3 && < 0.6.2+ other-modules: AbifTests+ , FastaTests+ , Generators+ ghc-options: -threaded -rtsopts -with-rtsopts=-N -Wall -Wincomplete-uni-patterns -Wincomplete-record-updates -Wimplicit-prelude -Wcompat -Wredundant-constraints -Wnoncanonical-monad-instances+ default-language: Haskell2010++source-repository head+ type: git+ location: https://github.com/hyraxbio/hyraxAbif
+ src/Examples.hs view
@@ -0,0 +1,10 @@+{-# LANGUAGE NoImplicitPrelude #-}+-- | Simple examples showing how to use this library+-- Also see "Hyrax.Abif.Generate" to see how a minimal AB1, that recall and phred can read, is generated+--+-- * "Examples.ReadAb1"+-- * "Examples.AddComment"+-- * "Examples.RemoveComments"+--+module Examples where+
+ src/Examples/AddComment.hs view
@@ -0,0 +1,29 @@+{-# LANGUAGE NoImplicitPrelude #-}+{-# LANGUAGE OverloadedStrings #-}++{-|+Description : Example showing how to add a comment to an existing AB1 file+Copyright : (c) HyraxBio, 2018+License : BSD3+Maintainer : andre@hyraxbio.co.za, andre@andrevdm.com++Example showing how to add a comment to an existing AB1 file.+See other examples in "Examples"+-}+module Examples.AddComment where++import Protolude++import qualified Hyrax.Abif.Read as H+import qualified Hyrax.Abif.Write as H++-- | Add a comment to an existing AB1 file+addComment :: IO ()+addComment = do+ abif' <- H.readAbif "example.ab1"++ case abif' of+ Left e -> putStrLn $ "error reading ABIF: " <> e+ Right abif -> do+ let modified = H.addDirectory abif $ H.mkComment "new comment"+ H.writeAbif "example.modified.ab1" modified
+ src/Examples/ReadAb1.hs view
@@ -0,0 +1,28 @@+{-# LANGUAGE NoImplicitPrelude #-}+{-# LANGUAGE OverloadedStrings #-}++{-|+Description : Example of reading a ABIF file+Copyright : (c) HyraxBio, 2018+License : BSD3+Maintainer : andre@hyraxbio.co.za, andre@andrevdm.com++Example of reading a ABIF file.+See other examples in "Examples"+-}+module Examples.ReadAb1 where++import Protolude++import qualified Hyrax.Abif.Read as H++-- | Read and print a ABIF file+addComment :: IO ()+addComment = do+ abif' <- H.readAbif "example.ab1"++ case abif' of+ Left e -> putStrLn $ "error reading ABIF: " <> e+ Right abif ->+ -- Print after removing the data, to make it readable+ print $ H.clearAbif abif
+ src/Examples/RemoveComments.hs view
@@ -0,0 +1,34 @@+{-# LANGUAGE NoImplicitPrelude #-}+{-# LANGUAGE OverloadedStrings #-}++{-|+Description : Example showing how to remove all comments from an ABIF file+Copyright : (c) HyraxBio, 2018+License : BSD3+Maintainer : andre@hyraxbio.co.za, andre@andrevdm.com++Example showing how to remove all comments from an ABIF file.+See other examples in "Examples"+-}+module Examples.RemoveComments where++import Protolude++import qualified Hyrax.Abif as H+import qualified Hyrax.Abif.Read as H+import qualified Hyrax.Abif.Write as H++-- | Remove all comments from an existing file+removeComments :: IO ()+removeComments = do+ abif' <- H.readAbif "example.ab1"++ case abif' of+ Left e -> putStrLn $ "error reading ABIF: " <> e+ Right abif -> do+ let modified = abif { H.aDirs = filter noComments $ H.aDirs abif }+ H.writeAbif "example.modified.ab1" modified++ where+ noComments :: H.Directory -> Bool+ noComments dir = H.dTagName dir /= "CMNT"
+ src/Hyrax/Abif.hs view
@@ -0,0 +1,124 @@+{-# LANGUAGE NoImplicitPrelude #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE MultiWayIf #-}++{-|+Module : Hyax.Abif+Description : Core AB1 types +Copyright : (c) HyraxBio, 2018+License : BSD3+Maintainer : andre@hyraxbio.co.za, andre@andrevdm.com+Stability : beta++This module contains the core types for working with AB1 files.++See++ * <https://github.com/hyraxbio/hyraxAbif/#readme Source code on github>++ * <http://www6.appliedbiosystems.com/support/software_community/ABIF_File_Format.pdf The ABIF spec>+-}+module Hyrax.Abif+ ( Abif (..)+ , Header (..)+ , Directory (..)+ , ElemType (..)+ , getElemType+ , describeElemType+ ) where++import Protolude+import qualified Data.ByteString.Lazy as BSL+++-- | A single ABIF+data Abif = Abif { aHeader :: !Header+ , aRootDir :: !Directory+ , aDirs :: ![Directory]+ } deriving (Show, Eq)+++-- | ABIF header+data Header = Header { hName :: !Text+ , hVersion :: !Int+ } deriving (Show, Eq)++-- | ABIF directory entry.+-- The 'dData' field contains the data for the entry+data Directory = Directory { dTagName :: !Text -- ^ Tag name+ , dTagNum :: !Int -- ^ Tag number, see e.g. how DATA entries use this+ , dElemType :: !ElemType -- ^ Type of an element+ , dElemTypeCode :: !Int -- ^ Integer value of 'dElemType'+ , dElemTypeDesc :: !Text -- ^ Description of 'dElemType'+ , dElemSize :: !Int -- ^ Size in bytes of each element+ , dElemNum :: !Int -- ^ Number of elements in the data. See the spec per data type. E.g. for a string this is the number of characters+ , dDataSize :: !Int -- ^ Number of bytes in the data+ , dDataOffset :: !Int -- ^ Offset of this directory entry's data in the file. For data that is four+ -- bytes or less, the data itself is stored in this field.+ -- This value will be recalculated when writing an ABIF so you do not need to manually set it.+ , dData :: !BSL.ByteString -- ^ The entry's data+ , dDataDebug :: ![Text] -- ^ Optinal debug data, populated by 'Hyrax.Abif.Read.getDebug' when a ABIF is parsed+ } deriving (Show, Eq)+++-- | Type of the elements in a directory entry. See the spec for details on each type if required.+data ElemType = ElemUnknown+ | ElemCustom+ | ElemByte+ | ElemChar+ | ElemWord+ | ElemShort+ | ElemLong+ | ElemFloat+ | ElemDouble+ | ElemDate+ | ElemTime+ | ElemPString+ | ElemCString+ | ElemThumb+ | ElemBool+ | ElemRationalUnsupported+ | ElemBCDUnsupported+ | ElemPointUnsupported+ | ElemRectUnsupported+ | ElemVPointUnsupported+ | ElemVRectUnsupported+ | ElemTagUnsupported+ | ElemDeltaCompUnsupported+ | ElemLZWCompUnsupported+ | ElemCompressedDataUnsupported+ | ElemRoot+ deriving (Show, Eq)+++-- | Get an 'ElemType' from a elem type code+getElemType :: Int -> ElemType+getElemType e = fst $ describeElemType e++-- | Get the description for an 'ElemType'+describeElemType :: Int -> (ElemType, Text)+describeElemType 1 = (ElemByte, "byte")+describeElemType 2 = (ElemChar, "char")+describeElemType 3 = (ElemWord, "word")+describeElemType 4 = (ElemShort, "short")+describeElemType 5 = (ElemLong, "long")+describeElemType 7 = (ElemFloat, "float")+describeElemType 8 = (ElemDouble, "double")+describeElemType 10 = (ElemDate, "date")+describeElemType 11 = (ElemTime, "time")+describeElemType 18 = (ElemPString, "pString")+describeElemType 19 = (ElemCString, "cString")+describeElemType 12 = (ElemThumb, "thumb")+describeElemType 13 = (ElemBool, "bool")+describeElemType 6 = (ElemRationalUnsupported, "rational (*unsupported*)")+describeElemType 9 = (ElemBCDUnsupported, "BCD (*unsupported*)")+describeElemType 14 = (ElemPointUnsupported, "point (*unsupported*)")+describeElemType 15 = (ElemRectUnsupported, "rect (*unsupported*)")+describeElemType 16 = (ElemVPointUnsupported, "vPoint (*unsupported*)")+describeElemType 17 = (ElemVRectUnsupported, "vRect (*unsupported*)")+describeElemType 20 = (ElemTagUnsupported, "Tag (*unsupported*)")+describeElemType 128 = (ElemDeltaCompUnsupported, "deltaComp (*unsupported*)")+describeElemType 256 = (ElemLZWCompUnsupported, "LZWComp (*unsupported*)")+describeElemType 384 = (ElemCompressedDataUnsupported, "Compressed Data (*unsupported*)")+describeElemType 1023 = (ElemRoot, "root")+describeElemType v = if v >= 1024 then (ElemCustom, "custom") else (ElemUnknown, "unknown")
+ src/Hyrax/Abif/Fasta.hs view
@@ -0,0 +1,50 @@+{-# LANGUAGE NoImplicitPrelude #-}+{-# LANGUAGE OverloadedStrings #-}++{-|+Module : Hyax.Abif.Fasta+Description : Read a FASTA file+Copyright : (c) HyraxBio, 2018+License : BSD3+Maintainer : andre@hyraxbio.co.za, andre@andrevdm.com+Stability : beta++Functionality for reading FASTA files+-}+module Hyrax.Abif.Fasta+ ( Fasta (..)+ , parseFasta+ ) where++import Protolude+import qualified Data.Text as Txt++-- | FASTA data+data Fasta = Fasta { fastaName :: !Text -- ^ Name+ , fastaRead :: !Text -- ^ Data+ } deriving (Show, Eq)+++-- | Parse the data for a single FASTA into a list of 'Fasta' values.+-- Single and multi-line FASTAs are supported.+-- Used by "Hyrax.Abif.Generate" to read weighted-FASTAs+parseFasta :: Text -> Either Text [Fasta]+parseFasta s =+ reverse <$> go (Txt.lines s) Nothing "" []++ where+ go :: [Text] -> Maybe Text -> Text -> [Fasta] -> Either Text [Fasta]+ go (line:lines) (Just name) read acc =+ if Txt.take 1 line /= ">"+ then go lines (Just name) (read <> line) acc+ else go lines (Just $ Txt.drop 1 line) "" (Fasta (Txt.strip name) read : acc)+ go (line:lines) Nothing _read acc =+ if Txt.take 1 line == ">"+ then go lines (Just $ Txt.strip . Txt.drop 1 $ line) "" acc+ else Left "Expecting name"+ go [] Nothing _ acc =+ Right acc+ go [] (Just _name) "" _acc =+ Left "Expecting read"+ go [] (Just name) read acc =+ Right $ Fasta (Txt.strip name) read : acc
+ src/Hyrax/Abif/Generate.hs view
@@ -0,0 +1,320 @@+{-# LANGUAGE NoImplicitPrelude #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE TupleSections #-}+{-|+Module : Hyax.Abif.Generate+Description : Generate AB1 from a weighted FASTA+Copyright : (c) HyraxBio, 2018+License : BSD3+Maintainer : andre@hyraxbio.co.za, andre@andrevdm.com++Functionality for generating AB1 files from an input FASTA. These AB1s are supported by both PHRED and recall,+if you are using other software you may need to add additional required sections.++= Weighted reads++The input FASTA files have "weighted" reads. The name for each read is an value between 0 and 1+ which specifies the height of the peak relative to a full peak. +++== Single read++The most simple example is a single FASTA with a single read with a weight of 1++@+> 1+ACTG+@++<<docs/eg_actg.png>>++The chromatogram for this AB1 shows perfect traces for the input `ACTG` nucleotides with a full height peak.+++== Mixes & multiple reads ++The source FASTA can have multiple reads, which results in a chromatogram with mixes++@+> 1+ACAG+> 0.3+ACTG+@++<<docs/eg_acag_acgt_mix.png>>++There is an `AT` mix at the third nucleotide. The first read has a weight of 1 and the second a weight of 0.3.+The chromatogram shows the mix and the `T` with a lower peak (30% of the `A` peak)++== Summing weights++ - The weigh of a read specifies the intensity of the peak from 0 to 1. + - Weights for each position are added to a maximum of 1 per nucleotide+ - You can use `_` as a "blank" nucleotide, in which only the nucleotides from other reads will be considered++E.g.++@+> 1+ACAG+> 0.3+_GT+> 0.2+_G+@++<<docs/eg_multi_mix.png>>++See README.md for additional details and examples+-}+module Hyrax.Abif.Generate+ ( generateAb1s+ , generateAb1+ , readWeightedFasta+ , iupac+ , unIupac+ ) where++import Protolude+import qualified Data.Text as Txt+import qualified Data.Text.Encoding as TxtE+import qualified Data.List as Lst+import qualified Data.Binary.Put as B+import qualified Data.ByteString as BS+import qualified Data.ByteString.Lazy as BSL+import qualified System.FilePath as FP+import System.FilePath ((</>))+import qualified System.Directory as Dir++import Hyrax.Abif+import Hyrax.Abif.Write+import Hyrax.Abif.Fasta++data TraceData = TraceData { trData09G :: ![Int16]+ , trData10A :: ![Int16]+ , trData11T :: ![Int16]+ , trData12C :: ![Int16]+ , trValsPerBase :: !Int+ , trFasta :: !Text+ } deriving (Show)++-- | Generate a set of AB1s. One for every FASTA found in the source directory+generateAb1s :: FilePath -> FilePath -> IO ()+generateAb1s source dest = do+ Dir.createDirectoryIfMissing True dest+ weighted <- readWeightedFastas source++ case weighted of+ Left e -> putText e+ Right rs -> do+ let ab1s = (\(n, r) -> (n, generateAb1 (n, r))) <$> rs+ traverse_ (\(name, ab1) -> BS.writeFile (dest </> Txt.unpack name <> ".ab1") $ BSL.toStrict ab1) ab1s+++-- | Create the 'ByteString' data for an AB1 given the data from a weighted FASTA (see 'readWeightedFasta')+generateAb1 :: (Text, [(Double, Text)]) -> BSL.ByteString+generateAb1 (fName, sourceFasta) = + let+ tr = generateTraceData sourceFasta+ valsPerBase = trValsPerBase tr+ generatedFastaLen = (Txt.length $ trFasta tr)++ -- The point that is the peak of the trace, i.e. mid point of trace for a single base+ midPeek = valsPerBase `div` 2+ -- Get the peak locations for all bases+ peakLocations = take generatedFastaLen [midPeek, valsPerBase + midPeek..]++ -- Sample name (from the FASTA name)+ sampleName = fst . Txt.breakOn "_" $ fName++ -- Create the ABIF directories+ dirs = [ mkData 9 $ trData09G tr -- G+ , mkData 10 $ trData10A tr -- A+ , mkData 11 $ trData11T tr -- T+ , mkData 12 $ trData12C tr -- C+ , mkBaseOrder BaseG BaseA BaseT BaseC -- Base order, should be GATC for 3500+ , mkLane 1 -- Lane or capliary number+ , mkCalledBases $ trFasta tr -- Called bases+ , mkMobilityFileName 1 "KB_3500_POP7_BDTv3.mob" -- Mobility file name+ , mkMobilityFileName 2 "KB_3500_POP7_BDTv3.mob" -- Mobility file name+ , mkPeakLocations $ fromIntegral <$> peakLocations -- Peak locations+ , mkDyeSignalStrength 53 75 79 48 -- Signal strength per dye+ , mkSampleName sampleName -- Sample name+ , mkComment "Generated by HyraxBio AB1 generator"+ ]++ -- The ABIF+ abif = Abif { aHeader = mkHeader+ , aRootDir = mkRoot+ , aDirs = dirs+ }+ + in+ -- Generate the data+ B.runPut (putAbif abif)+++-- | Generate the traces for the AB1 from the parsed weighted FASTA+generateTraceData :: [(Double, Text)] -> TraceData+generateTraceData weighted =+ let+ weightedNucs' = (\(w, ns) -> (w,) . unIupac <$> Txt.unpack ns) <$> weighted+ weightedNucs = Lst.transpose weightedNucs'+ + -- Values for a base that was present. This defines the shape of the chromatogram curve, and defines the number of values per base+ curve = [0, 0, 128, 512, 1024, 1024, 512, 128, 0, 0]+ valsPerBase = length curve++ -- Create the G, A, T and C traces+ data09G = concat $ getWeightedTrace curve 'G' <$> weightedNucs+ data10A = concat $ getWeightedTrace curve 'A' <$> weightedNucs+ data11T = concat $ getWeightedTrace curve 'T' <$> weightedNucs+ data12C = concat $ getWeightedTrace curve 'C' <$> weightedNucs++ -- Create fasta sequence for the trace+ fastaSeq = concat <$> (snd <<$>> weightedNucs)+ fasta = Txt.pack $ iupac fastaSeq+ in + TraceData { trData09G = data09G+ , trData10A = data10A+ , trData11T = data11T+ , trData12C = data12C+ , trFasta = fasta+ , trValsPerBase = valsPerBase+ }++ where+ getWeightedTrace :: [Int] -> Char -> [(Double, [Char])] -> [Int16]+ getWeightedTrace curve nuc ws =+ let+ found = filter ((nuc `elem`) . snd) ws+ score' = foldl' (+) 0 $ fst <$> found+ score = min 1 . max 0 $ score'+ wave = floor . (score *) . fromIntegral <$> curve+ in+ wave+++-- | Read a weighted FASTA file. See the module comments for the expected format.+-- See the module documentation for details on the format of the weighted FASTA +--+-- e.g. weighted FASTA+--+-- @+-- > 1+-- ACAG+-- > 0.3+-- _GT+-- > 0.2+-- _G+-- @+--+--+-- The result data has the type+-- +-- @+-- ('Text', [('Double', 'Text')])+-- ^ ^ ^+-- | | |+-- file name -------------+ | +---- read +-- | +-- +---- weight+-- @+--+readWeightedFasta :: ByteString -> Either Text [(Double, Text)]+readWeightedFasta fastaData = + case parseFasta $ TxtE.decodeUtf8 fastaData of+ Left e -> Left e+ Right fs -> getWeightedFasta fs++ where+ getWeightedFasta :: [Fasta] -> Either Text [(Double, Text)]+ getWeightedFasta fs = + case sequenceA $ readWeighted <$> fs of+ Left e -> Left e+ Right r -> Right r++ readWeighted :: Fasta -> Either Text (Double, Text)+ readWeighted (Fasta hdr dta) =+ case (readMaybe . Txt.unpack $ hdr :: Maybe Double) of+ Just weight -> Right (min 1 . max 0 $ weight, Txt.strip dta)+ Nothing -> Left $ "Invalid header reading, expecting numeric weight, got: " <> hdr++ ++-- | Read all FASTA files in a directory+readWeightedFastas :: FilePath -> IO (Either Text [(Text, [(Double, Text)])])+readWeightedFastas source = do+ files <- filter (Txt.isSuffixOf ".fasta" . Txt.pack) <$> getFiles source+ let names = Txt.pack . FP.takeBaseName <$> files+ contents <- traverse BS.readFile files+ + case sequenceA $ readWeightedFasta <$> contents of+ Left e -> pure . Left $ e+ Right rs -> pure . Right $ zip names rs++ +-- | Find all files in a directory+getFiles :: FilePath -> IO [FilePath]+getFiles p = do+ entries <- (p </>) <<$>> Dir.listDirectory p+ filterM Dir.doesFileExist entries+++-- | Convert a IUPAC ambiguity code to the set of nucleotides it represents+unIupac :: Char -> [Char]+unIupac c =+ case c of+ 'T' -> "T"+ 'C' -> "C"+ 'A' -> "A"+ 'G' -> "G"+ + 'U' -> "T"+ 'M' -> "AC"+ 'R' -> "AG"+ 'W' -> "AT"+ 'S' -> "CG"+ 'Y' -> "CT"+ 'K' -> "GT"+ 'V' -> "ACG"+ 'H' -> "ACT"+ 'D' -> "AGT"+ 'B' -> "CGT"+ 'N' -> "GATC"+ + 'X' -> "GATC"+ _ -> ""+++-- | Given a set of nucleotides get the IUPAC ambiguity code+iupac :: [[Char]] -> [Char]+iupac ns =+ go <$> ns++ where+ go cs =+ let+ a = 'A' `elem` cs+ c = 'C' `elem` cs+ g = 'G' `elem` cs+ t = 'T' `elem` cs+ in+ case (a, c, g, t) of+ (True, False, False, False) -> 'A'+ (False, True, False, False) -> 'C'+ (False, False, True, False) -> 'G'+ (False, False, False, True ) -> 'T'+ (True, True, False, False) -> 'M'+ (True, False, True, False) -> 'R'+ (True, False, False, True ) -> 'W'+ (False, True, True, False) -> 'S'+ (False, True, False, True ) -> 'Y'+ (False, False, True, True ) -> 'K'+ (True, True, True, False) -> 'V'+ (True, True, False, True ) -> 'H'+ (True, False, True, True ) -> 'D'+ (False, True, True, True ) -> 'B'+ (True, True, True, True ) -> 'N'+ _ -> '_'
+ src/Hyrax/Abif/Read.hs view
@@ -0,0 +1,249 @@+{-# LANGUAGE NoImplicitPrelude #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE MultiWayIf #-}++{-|+Module : Hyax.Abif.Read+Description : Read and parse AB1 files+Copyright : (c) HyraxBio, 2018+License : BSD3+Maintainer : andre@hyraxbio.co.za, andre@andrevdm.com+Stability : beta++Functionality for reading and parsing AB1 files++e.g.++@+abif' <- readAbif "example.ab1"++case abif' of+ Left e -> putStrLn $ "error reading ABIF: " <> e+ Right abif -> print $ clearAbif abif+@+-}+module Hyrax.Abif.Read+ ( readAbif+ , getAbif+ , clear+ , clearAbif+ , getDebug+ , getPString+ , getCString+ , getHeader+ , getRoot+ , getDirectories+ , getDirectory+ ) where++import Protolude+import qualified Data.Text as Txt+import qualified Data.Text.Encoding as TxtE+import qualified Data.Binary as B+import qualified Data.Binary.Get as B+import qualified Data.ByteString.Lazy as BSL+import Control.Monad.Fail (fail)++import Hyrax.Abif+++-- | Read and parse an AB1 file+readAbif :: FilePath -> IO (Either Text Abif)+readAbif path = getAbif <$> BSL.readFile path+++-- | Parse an AB1 from a 'ByteString'+getAbif :: BSL.ByteString -> Either Text Abif+getAbif bs = do+ (header, rootDir) <- case B.runGetOrFail (getRoot bs) bs of+ Right (_, _, x) -> pure x+ Left (_, _, e) -> Left ("Error reading root: " <> Txt.pack e)++ let dirBytes = BSL.drop (fromIntegral $ dDataOffset rootDir) bs+ + ds <- case B.runGetOrFail (getDirectories bs [] $ dElemNum rootDir) dirBytes of+ Right (_, _, x) -> pure x+ Left (_, _, e) -> Left ("Error reading " <> show (dElemNum rootDir) <> " directories (at " <> show (dDataOffset rootDir) <> "): " <> Txt.pack e)+ + pure $ Abif header rootDir ds+++-- | Removes all data from the ABIF's directories+clearAbif :: Abif -> Abif+clearAbif a = a { aRootDir = clear $ aRootDir a+ , aDirs = clear <$> aDirs a+ }+++-- | Removes all data from a directory entry. This will probably only be useful when trying to show an ABIF value+clear :: Directory -> Directory+clear d = d { dData = "" }+++-- | Populate the directory entry with debug data (into 'dDataDebug').+-- This is done for selected types only, e.g. for strings so that printing the structure will display+-- readable/meaningfull info+getDebug :: Directory -> Directory+getDebug d =+ let bsAtOffset = dData d in+ + case dElemType d of+ -- Strings have a count = number of chars, not number of "strings"+ ElemPString ->+ if dDataSize d <= 4+ then d { dDataDebug = [TxtE.decodeUtf8 . BSL.toStrict . BSL.drop 1 . BSL.take (fromIntegral $ dDataSize d) $ dData d] }+ else d { dDataDebug = [B.runGet (lbl getPString) bsAtOffset] }++ -- Strings have a count = number of chars, not number of "strings"+ ElemCString ->+ if dDataSize d <= 4+ then d { dDataDebug = [TxtE.decodeUtf8 . BSL.toStrict . BSL.take (fromIntegral $ dDataSize d - 1) $ dData d] }+ else d { dDataDebug = [B.runGet (lbl . getCString $ dDataSize d) bsAtOffset] }++ y ->+ -- For non-array entries+ if dElemNum d == 1+ then + case y of+ ElemDate -> + flip B.runGet (dData d) $ lbl $ do+ yy <- B.getInt16be+ mt <- B.getInt8+ dt <- B.getInt8+ pure d { dDataDebug = [show yy <> "/" <> show mt <> "/" <> show dt]}+ + ElemTime ->+ flip B.runGet (dData d) $ lbl $ do+ hr <- B.getInt8+ mn <- B.getInt8+ sc <- B.getInt8+ ss <- B.getInt8+ pure $ d { dDataDebug = [show hr <> ":" <> show mn <> ":" <> show sc <> "." <> show ss] }+ + ElemLong ->+ flip B.runGet (dData d) $ lbl $ do+ x <- B.getInt32be+ pure $ d { dDataDebug = [show x] }+ + ElemShort ->+ flip B.runGet (dData d) $ lbl $ do+ x <- B.getInt16be+ pure $ d { dDataDebug = [show x] }+ + ElemFloat ->+ flip B.runGet (dData d) $ lbl $ do+ x <- B.getFloatbe+ pure $ d { dDataDebug = [show x] }+ + ElemWord ->+ flip B.runGet (dData d) $ lbl $ do+ x <- B.getInt8+ pure $ d { dDataDebug = [show x] }+ + ElemChar ->+ flip B.runGet (dData d) $ lbl $ do+ x <- B.getWord8+ let c = BSL.pack [x]+ pure $ d { dDataDebug = [TxtE.decodeUtf8 . BSL.toStrict $ c] }+ + _ -> d+ else+ case y of+ ElemChar -> -- Array of chars can be treated as a string+ flip B.runGet (dData d) $ lbl $ do+ cs <- readArray B.getWord8+ let c = BSL.pack cs+ pure $ d { dDataDebug = [TxtE.decodeUtf8 . BSL.toStrict $ c] }++ --ElemShort ->+ -- flip B.runGet (dData d) $ lbl $ do+ -- xs <- readArray B.getInt16be+ -- pure $ d { dDataDebug = [show xs] }++ _ -> d -- Do nothing++ where+ lbl = B.label $ "Reading " <> show (dElemTypeDesc d) <> " data size=" <> show (dDataSize d) <> " dir entry=" <> Txt.unpack (dTagName d) <> " cached data size=" <> (show . BSL.length $ dData d) <> ". "++ readArray :: B.Get n -> B.Get [n]+ readArray getFn = do+ e <- B.isEmpty+ if e then return []+ else do+ c <- getFn+ cs <- readArray getFn+ pure (c:cs)+++-- | Parse a 'ElemPString'+getPString :: B.Get Text+getPString = do+ sz <- fromIntegral <$> B.getInt8+ TxtE.decodeUtf8 <$> B.label ("PString length=" <> show sz <> ".") (B.getByteString sz)+++-- | Parse a 'ElemCString'+getCString :: Int -> B.Get Text+getCString sz = + TxtE.decodeUtf8 <$> B.getByteString (sz - 1)+++-- | Parse the ABIF 'Header'+getHeader :: B.Get Header+getHeader = + Header <$> (TxtE.decodeUtf8 <$> B.getByteString 4)+ <*> (fromIntegral <$> B.getInt16be)+++-- | Parse the root ('Header' and 'Directory')+getRoot :: BSL.ByteString -> B.Get (Header, Directory)+getRoot bs = do+ h <- getHeader+ rd <- getDirectory bs+ pure (h, rd)+++-- | Parse a single 'Directory' entry and read its data+getDirectory :: BSL.ByteString -> B.Get Directory+getDirectory bs = do+ tagName <- TxtE.decodeUtf8 <$> B.getByteString 4+ tagNum <- fromIntegral <$> B.getInt32be+ typeCode <- fromIntegral <$> B.getInt16be+ elemSize <- fromIntegral <$> B.getInt16be+ elemNum <- fromIntegral <$> B.getInt32be+ dataSize <- fromIntegral <$> B.getInt32be+ offsetDataBytes <- B.lookAhead $ B.getLazyByteString 4+ dataOffset <- fromIntegral <$> B.getInt32be++ -- Read the data+ -- Data that is 4 bytes or less is stored in the offset field+ dataBytes <- if dataSize <= 4+ then pure $ BSL.take (fromIntegral dataSize) offsetDataBytes+ else case B.runGetOrFail (B.getLazyByteString $ fromIntegral dataSize) $ BSL.drop (fromIntegral dataOffset) bs of+ Right (_, _, x) -> pure x+ Left (_, _, e) -> fail $ "error reading data (" <> show dataSize <> " bytes starting at " <> show dataOffset <> ") for directory entry '" <> Txt.unpack tagName <> "': " <> e++ let (elemType, elemCode) = describeElemType typeCode+ pure Directory { dTagName = tagName + , dTagNum = tagNum + , dElemTypeCode = typeCode + , dElemTypeDesc = elemCode + , dElemType = elemType + , dElemSize = elemSize + , dElemNum = elemNum + , dDataSize = dataSize + , dDataOffset = dataOffset + , dData = dataBytes + , dDataDebug = []+ } +++-- | Parse all the directoy entries+getDirectories :: BSL.ByteString -> [Directory] -> Int -> B.Get [Directory]+getDirectories _ acc 0 = pure acc+getDirectories bs acc more = do+ d <- getDirectory bs+ B.skip 4 -- Skip the reserved field+ getDirectories bs (acc <> [d]) (more - 1)++
+ src/Hyrax/Abif/Write.hs view
@@ -0,0 +1,334 @@+{-# LANGUAGE NoImplicitPrelude #-}+{-# LANGUAGE OverloadedStrings #-}++{-|+Module : Hyax.Abif.Write+Description : Functionality for writing AB1 files+Copyright : (c) HyraxBio, 2018+License : BSD3+Maintainer : andre@hyraxbio.co.za, andre@andrevdm.com+Stability : beta++Functionality for writing AB1 files.+See 'Hyrax.Abif.Generate.generateAb1' for an example of how to create an 'Ab1'+-}+module Hyrax.Abif.Write+ ( createAbifBytes+ , writeAbif+ , putAbif+ , putTextStr+ , putHeader+ , putDirectory+ , mkHeader+ , mkRoot+ , mkData+ , mkComment+ , mkSampleName+ , mkBaseOrder+ , mkLane+ , mkCalledBases+ , mkMobilityFileName+ , mkDyeSignalStrength+ , mkPeakLocations+ , addDirectory+ , Base (..)+ ) where++import Protolude+import qualified Data.Text as Txt+import qualified Data.Text.Encoding as TxtE+import qualified Data.Binary as B+import qualified Data.Binary.Put as B+import qualified Data.ByteString as BS+import qualified Data.ByteString.Lazy as BSL++import Hyrax.Abif++-- | Used to specify the base order for the FWO directry entry, see 'mkBaseOrder'+data Base = BaseA | BaseC | BaseG | BaseT+++-- | Write an 'Abif' to a 'ByteString'+createAbifBytes :: Abif -> BSL.ByteString+createAbifBytes ab1 =+ B.runPut (putAbif ab1)++ +-- | Write an 'Abif' to a file+writeAbif :: FilePath -> Abif -> IO ()+writeAbif destPath ab1 = do+ let b = createAbifBytes ab1+ BS.writeFile destPath $ BSL.toStrict b++ +-- | Create the 'Abif' using "Data.Binary"+putAbif :: Abif -> B.Put+putAbif (Abif header root dirs) = do+ -- Data starts at offset 128+ let startDataOffset = 128+ -- Total data size+ let dataSize = foldl' (\acc i -> if i > 4 then acc + i else acc) 0 $ dDataSize <$> dirs+ + -- Write the header+ putHeader header++ -- Write the root directory entry+ putDirectory (startDataOffset + dataSize) $ root { dDataSize = 28 * length dirs+ , dElemNum = length dirs+ }++ -- Write 47 zero Int16 values as required by the spec+ traverse_ B.putInt16be $ replicate 47 0+ -- Write the data, for all data larger than four bytes. Data four bytes or less is stored+ -- in the offset field+ traverse_ (B.putLazyByteString . dData) $ filter (\d -> dDataSize d > 4) dirs+ -- Write the directory entries. + foldM_ writeDir startDataOffset dirs++ where+ writeDir offset dir = do+ putDirectory offset dir+ pure $ if dDataSize dir > 4+ then offset + dDataSize dir+ else offset+++-- | Write 'Text'+putTextStr :: Text -> B.Put+putTextStr t = B.putByteString $ TxtE.encodeUtf8 t+++-- | Write a 'ElemPString'+putPStr :: Text -> B.Put+putPStr t = do+ B.putInt8 . fromIntegral $ Txt.length t+ B.putByteString $ TxtE.encodeUtf8 t+++-- | Write a 'Header'+putHeader :: Header -> B.Put+putHeader h = do+ putTextStr $ hName h+ B.putInt16be . fromIntegral $ hVersion h+++-- | Write a 'Directory'+putDirectory :: Int -> Directory -> B.Put+putDirectory dirOffset d = do+ let name = Txt.justifyLeft 4 ' ' . Txt.take 4 $ dTagName d+ putTextStr name+ B.putInt32be . fromIntegral $ dTagNum d+ B.putInt16be . fromIntegral $ dElemTypeCode d+ B.putInt16be . fromIntegral $ dElemSize d+ B.putInt32be . fromIntegral $ dElemNum d+ B.putInt32be . fromIntegral $ dDataSize d++ -- data with a size >= 4 are written in the offset+ if dDataSize d > 4+ then B.putInt32be . fromIntegral $ dirOffset+ else B.putLazyByteString . BSL.take 4 $ dData d <> "\0\0\0\0"++ B.putInt32be 0 -- reserved / datahandle+++-- | Create a 'Header'+mkHeader :: Header+mkHeader =+ Header { hName = "ABIF"+ , hVersion = 101+ }+++-- | Create the root 'Directory' entry+mkRoot :: Directory+mkRoot = + Directory { dTagName = "tdir"+ , dTagNum = 1+ , dElemTypeCode = 1023+ , dElemTypeDesc = "root"+ , dElemType = ElemRoot+ , dElemSize = 28+ , dDataOffset = 0+ , dDataDebug = []+ , dData = ""+ , dDataSize = 0+ , dElemNum = 0+ }+ ++-- | Create a comment 'Directory' entry and 'ElemPString' data+mkComment :: Text -> Directory+mkComment comment' = + let comment = B.runPut . putPStr $ comment' in++ Directory { dTagName = "CMNT" -- Comment+ , dTagNum = 1+ , dElemTypeCode = 18+ , dElemTypeDesc = "pString"+ , dElemType = ElemPString+ , dElemSize = 1+ , dElemNum = 1+ , dDataOffset = 0+ , dDataDebug = []+ , dData = comment+ , dDataSize = fromIntegral (BSL.length comment)+ } +++-- | Create a sample name (SMPL) 'Directory' entry and 'ElemPString' data+mkSampleName :: Text -> Directory+mkSampleName sampleName' =+ let sampleName = B.runPut . putPStr $ sampleName' in+ Directory { dTagName = "SMPL" -- Sample name+ , dTagNum = 1+ , dElemTypeCode = 18+ , dElemTypeDesc = "pString"+ , dElemType = ElemPString+ , dElemSize = 1+ , dElemNum = 10+ , dDataOffset = 0+ , dDataDebug = []+ , dData = sampleName+ , dDataSize = fromIntegral (BSL.length sampleName)+ }++-- | Create a base order (FWO_) 'Directory' entry data+mkBaseOrder :: Base -> Base -> Base -> Base -> Directory+mkBaseOrder w x y z =+ Directory { dTagName = "FWO_" -- Base order+ , dTagNum = 1+ , dElemTypeCode = 2+ , dElemTypeDesc = "char"+ , dElemType = ElemChar+ , dElemSize = 1+ , dDataOffset = 0+ , dDataDebug = []+ , dData = getBase w <> getBase x <> getBase y <> getBase z+ , dDataSize = 4+ , dElemNum = 4+ }+ where+ getBase BaseA = "A"+ getBase BaseC = "C"+ getBase BaseG = "G"+ getBase BaseT = "T"+++-- | Create a lane (LANE) 'Directory' entry and data+mkLane :: Int16 -> Directory+mkLane lane =+ Directory { dTagName = "LANE" -- Lane or capliary number+ , dTagNum = 1+ , dElemTypeCode = 4+ , dElemTypeDesc = "short"+ , dElemType = ElemShort+ , dElemSize = 2+ , dElemNum = 1+ , dDataSize = 2+ , dDataOffset = 0+ , dData = B.runPut $ B.putInt16be lane+ , dDataDebug = []+ }+++-- | Create a called bases (PBAS) 'Directory' entry and data+mkCalledBases :: Text -> Directory+mkCalledBases fasta = + let+ generatedFastaLen = Txt.length fasta+ pbas = BSL.fromStrict . TxtE.encodeUtf8 $ fasta+ in+ Directory { dTagName = "PBAS" -- Called bases+ , dTagNum = 1+ , dElemTypeCode = 2+ , dElemTypeDesc = "char"+ , dElemType = ElemChar+ , dElemSize = 1+ , dDataOffset = 0+ , dDataDebug = []+ , dData = pbas+ , dDataSize = generatedFastaLen+ , dElemNum = generatedFastaLen + }+++-- | Create a mobility file name (PDMF) 'Directory' entry and 'ElemPString' data+mkMobilityFileName :: Int -> Text -> Directory+mkMobilityFileName tagNum fileName =+ let pdfm = B.runPut $ putPStr fileName in+ Directory { dTagName = "PDMF" -- Mobility file name+ , dTagNum = tagNum+ , dElemTypeCode = 18+ , dElemTypeDesc = "pString"+ , dElemType = ElemPString+ , dElemSize = 1+ , dDataOffset = 0+ , dDataDebug = []+ , dData = pdfm+ , dDataSize = fromIntegral (BSL.length pdfm)+ , dElemNum = fromIntegral (BSL.length pdfm) + }+++-- | Create a signal strength (S/N%) 'Directory' entry and data+mkDyeSignalStrength :: Int16 -> Int16 -> Int16 -> Int16 -> Directory+mkDyeSignalStrength w x y z =+ let sigStrength = B.runPut $ do+ B.putInt16be w+ B.putInt16be x+ B.putInt16be y+ B.putInt16be z+ in+ Directory { dTagName = "S/N%" -- Signal strength per dye+ , dTagNum = 1+ , dElemTypeCode = 4+ , dElemTypeDesc = "short"+ , dElemType = ElemShort+ , dElemSize = 2+ , dElemNum = 4+ , dDataOffset = 0+ , dDataDebug = []+ , dData = sigStrength+ , dDataSize = fromIntegral (BSL.length sigStrength)+ }+++-- | Create a peak locations (PLOC) 'Directory' entry and array of 'ElemShort' data+mkPeakLocations :: [Int16] -> Directory+mkPeakLocations locs =+ let peakLocations = B.runPut $ traverse_ B.putInt16be locs in+ Directory { dTagName = "PLOC" -- Peak locations+ , dTagNum = 1+ , dElemTypeCode = 4+ , dElemTypeDesc = "short"+ , dElemType = ElemShort+ , dElemSize = 2+ , dDataOffset = 0+ , dDataDebug = []+ , dData = peakLocations+ , dDataSize = fromIntegral $ BSL.length peakLocations+ , dElemNum = length locs+ }+++-- | Create a data (DATA) 'Directory' entry and array of 'ElemShort' data+mkData :: Int -> [Int16] -> Directory+mkData tagNum ds =+ let ds' = B.runPut $ traverse_ B.putInt16be ds in+ Directory { dTagName = "DATA"+ , dTagNum = tagNum+ , dElemTypeCode = 4+ , dElemTypeDesc = "short"+ , dElemType = ElemShort+ , dElemSize = 2+ , dDataOffset = 0+ , dDataDebug = []+ , dData = ds'+ , dDataSize = fromIntegral (BSL.length ds')+ , dElemNum = length ds+ }++-- | Add a directory to an 'Abif'+addDirectory :: Abif -> Directory -> Abif+addDirectory abif dir =+ abif { aDirs = aDirs abif <> [dir] }
+ test/AbifTests.hs view
@@ -0,0 +1,119 @@+{-# LANGUAGE TemplateHaskell #-}+{-# LANGUAGE NoImplicitPrelude #-}+{-# LANGUAGE OverloadedStrings #-}++module AbifTests (tests) where+++import Protolude+import qualified Data.List as Lst+import qualified Data.Text as Txt+import qualified Data.Text.Encoding as TxtE+import qualified Data.Binary as B+import qualified Data.Binary.Get as B+import qualified Data.ByteString.Lazy as BSL+import Hedgehog++import qualified Hyrax.Abif as H+import qualified Hyrax.Abif.Read as H+import qualified Hyrax.Abif.Write as H+import qualified Hyrax.Abif.Generate as H+import Generators++-- | Test that an ab1 (write, read, write, read) results in the original data+prop_roundtrip :: Property+prop_roundtrip = property $ do + fdata <- forAll $ genFastaData 2+ let fasta = TxtE.encodeUtf8 $ toFastaTxt False fdata++ wfasta <- evalEither $ H.readWeightedFasta fasta+ let ab1Written1 = H.generateAb1 ("test", wfasta)+ ab1Read1 <- evalEither $ H.getAbif ab1Written1++ let ab1Written2 = H.createAbifBytes ab1Read1+ ab1Read2 <- evalEither $ H.getAbif ab1Written2++ ab1Read1 === ab1Read2+ + +-- | Generate an ab1 from a fasta and then confirm that the generated peaks of the+-- chromatogram match the original fasta.+-- Note that we are only testing the simple/single fasta case (i.e. no mixes)+prop_readPeaks :: Property+prop_readPeaks = property $ do+ nucs <- nucsNoIupacGen+ let fasta = TxtE.encodeUtf8 $ "> 1\n" <> nucs++ wfasta <- evalEither $ H.readWeightedFasta fasta+ let ab1Written = H.generateAb1 ("test", wfasta)+ ab1 <- evalEither $ H.getAbif ab1Written++ -- Get the peak locations+ peaks <- evalEither $ readShorts <$> getDirEntry ab1 "PLOC" 1++ -- Get the val at the peak per channel+ chanG <- evalEither $ readDataPeaks ab1 9 peaks+ chanA <- evalEither $ readDataPeaks ab1 10 peaks+ chanT <- evalEither $ readDataPeaks ab1 11 peaks+ chanC <- evalEither $ readDataPeaks ab1 12 peaks++ -- Check that all the data is available+ length chanA === length peaks+ length chanC === length peaks+ length chanG === length peaks+ length chanT === length peaks++ -- Call the peaks+ let called = Txt.pack $ Lst.zipWith4 callPeaks chanA chanC chanG chanT++ -- Compare original fasta vs called+ called === nucs++ -- Compare to the sequence in the AB1+ pbas <- evalEither $ getDirEntry ab1 "PBAS" 1+ (TxtE.decodeUtf8 . BSL.toStrict $ pbas) === nucs++ where+ -- Highest peak wins, no iupac+ callPeaks :: Int -> Int -> Int -> Int -> Char+ callPeaks a c g t =+ case reverse $ Lst.sortOn snd [('A', a), ('C', c), ('G', g), ('T', t)] of+ ((n,_) : _) -> n+ _ -> '?'+ + readDataPeaks :: H.Abif -> Int -> [Int] -> Either Text [Int]+ readDataPeaks ab1 dirNum peaks =+ case readShorts <$> getDirEntry ab1 "DATA" dirNum of+ Left e -> Left e+ Right vs ->+ let valsAtPeaks = atMay vs <$> peaks in+ maybeToRight "peaks" $ sequenceA valsAtPeaks+ ++getDirEntry :: H.Abif -> Text -> Int -> Either Text BSL.ByteString+getDirEntry ab1 dirName dirNum =+ let r = filter (\d -> H.dTagNum d == dirNum && H.dTagName d == dirName) $ H.aDirs ab1 in+ case r of+ (a:_) -> Right $ H.dData a+ _ -> Left $ "No entry found for '" <> dirName <> "' " <> show dirNum+++readShorts :: BSL.ByteString -> [Int]+readShorts =+ let r = B.runGet $ readArray B.getWord16be in+ fromIntegral <<$>> r+++readArray :: B.Get n -> B.Get [n]+readArray getFn = do+ e <- B.isEmpty+ if e then return []+ else do+ c <- getFn+ cs <- readArray getFn+ pure (c:cs)+++tests :: IO Bool+tests =+ checkParallel $$discover
+ test/FastaTests.hs view
@@ -0,0 +1,27 @@+{-# LANGUAGE TemplateHaskell #-}+{-# LANGUAGE NoImplicitPrelude #-}+{-# LANGUAGE OverloadedStrings #-}++module FastaTests (tests) where++import Protolude+import Hedgehog+import qualified Hedgehog.Gen as Gen++import qualified Hyrax.Abif.Fasta as F+import Generators++-- | Test fasta parsing+prop_fasta :: Property+prop_fasta = property $ do+ withSpaces <- forAll $ Gen.element [True, False]+ fdata <- forAll $ genFastaData 8+ let ftxt = toFastaTxt withSpaces fdata+ let fasta = toFasta fdata+ let parsed = F.parseFasta ftxt+ parsed === Right fasta+ ++tests :: IO Bool+tests =+ checkParallel $$discover
+ test/Generators.hs view
@@ -0,0 +1,53 @@+{-# LANGUAGE NoImplicitPrelude #-}+{-# LANGUAGE OverloadedStrings #-}++module Generators where++import Protolude+import qualified Data.Text as Txt+import Hedgehog+import qualified Hedgehog.Gen as Gen+import qualified Hedgehog.Range as Range++import qualified Hyrax.Abif.Fasta as F++-- | Generate random set of nucleotides, including ambiguous ones+nucsGen :: GenT Identity Text+nucsGen = + Gen.text (Range.linear 1 1000) (Gen.element "ACGTMRWSYKVHDBNX")+++-- | Generate random set of nucleotides, no ambiguous ones+nucsNoIupacGen :: (Monad m) => PropertyT m Text+nucsNoIupacGen = + forAll $ Gen.text (Range.linear 1 1000) (Gen.element "ACGT")++ +-- | Generate data to construct a fasta "file"+-- List of tuples of (name, lines of reads)+-- maxReads controls the max number in the list+genFastaData :: Int -> Gen [(Int, [Text])]+genFastaData maxReads = + Gen.list (Range.linear 1 maxReads) $+ (,) <$> Gen.int (Range.linear 1 9999) <*> Gen.list (Range.linear 1 15) nucsGen+++-- | Convert fasta data to a single string+toFastaTxt :: Bool -> [(Int, [Text])] -> Text+toFastaTxt addSpaces is =+ Txt.intercalate "\n" $ toRead <$> is+ + where+ toRead (weight, ls) =+ let sp = if addSpaces then " " else "" in+ ">" <> sp <> show weight <> "\n" <> Txt.intercalate "\n" ls+++-- | Convert fasta data to the expected list of FASTA values+toFasta :: [(Int, [Text])] -> [F.Fasta]+toFasta is =+ toRead <$> is+ + where+ toRead (weight, ls) = F.Fasta (show weight) $ Txt.concat ls+
+ test/Tests.hs view
@@ -0,0 +1,21 @@+{-# LANGUAGE NoImplicitPrelude #-}++import Protolude+import qualified System.IO as IO+import System.Exit (exitFailure)++import qualified AbifTests+import qualified FastaTests++main :: IO ()+main = do+ IO.hSetBuffering IO.stdout IO.LineBuffering+ IO.hSetBuffering IO.stderr IO.LineBuffering++ results <- sequence [ AbifTests.tests+ , FastaTests.tests+ ]++ if and results+ then pass+ else exitFailure