diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -138,6 +138,14 @@ *Note that the reads do not need to be the same length.* + +### Reverse reads + +A weighted FASTA can represent a reverse read. To do this add a `R` suffix to the weight. +The data you enter should be entered as if it was a forward read. This data will be complemented +and reversed before writing to the ABIF + + --- #### Example FASTA - single file diff --git a/hyraxAbif.cabal b/hyraxAbif.cabal --- a/hyraxAbif.cabal +++ b/hyraxAbif.cabal @@ -1,6 +1,6 @@ cabal-version: 2.0 name: hyraxAbif -version: 0.2.3.5 +version: 0.2.3.7 synopsis: Modules for parsing, generating and manipulating AB1 files. homepage: https://github.com/hyraxbio/hyraxAbif/#readme license: BSD3 @@ -60,7 +60,7 @@ , 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 + , pretty-show >= 1.6.16 && < 1.9.0 , hscolour >= 1.24.4 && < 1.25.0 default-language: Haskell2010 diff --git a/src/Examples/ReadAb1.hs b/src/Examples/ReadAb1.hs --- a/src/Examples/ReadAb1.hs +++ b/src/Examples/ReadAb1.hs @@ -17,8 +17,8 @@ import qualified Hyrax.Abif.Read as H -- | Read and print a ABIF file -addComment :: IO () -addComment = do +readAbif :: IO () +readAbif = do abif' <- H.readAbif "example.ab1" case abif' of diff --git a/src/Hyrax/Abif.hs b/src/Hyrax/Abif.hs --- a/src/Hyrax/Abif.hs +++ b/src/Hyrax/Abif.hs @@ -18,6 +18,7 @@ * -} +{-! SECTION< abif_module !-} module Hyrax.Abif ( Abif (..) , Header (..) @@ -29,8 +30,10 @@ import Protolude import qualified Data.ByteString.Lazy as BSL +{-! SECTION> abif_module !-} +{-! SECTION< abif_type !-} -- | A single ABIF data Abif = Abif { aHeader :: !Header , aRootDir :: !Directory @@ -59,8 +62,10 @@ , 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) +{-! SECTION> abif_type !-} +{-! SECTION< abif_elemType !-} -- | Type of the elements in a directory entry. See the spec for details on each type if required. data ElemType = ElemUnknown | ElemCustom @@ -122,3 +127,4 @@ describeElemType 384 = (ElemCompressedDataUnsupported, "Compressed Data (*unsupported*)") describeElemType 1023 = (ElemRoot, "root") describeElemType v = if v >= 1024 then (ElemCustom, "custom") else (ElemUnknown, "unknown") +{-! SECTION> abif_elemType !-} diff --git a/src/Hyrax/Abif/Fasta.hs b/src/Hyrax/Abif/Fasta.hs --- a/src/Hyrax/Abif/Fasta.hs +++ b/src/Hyrax/Abif/Fasta.hs @@ -19,6 +19,7 @@ import Protolude import qualified Data.Text as Txt +{-! SECTION< fasta !-} -- | FASTA data data Fasta = Fasta { fastaName :: !Text -- ^ Name , fastaRead :: !Text -- ^ Data @@ -48,3 +49,4 @@ Left "Expecting read" go [] (Just name) read acc = Right $ Fasta (Txt.strip name) read : acc +{-! SECTION> fasta !-} diff --git a/src/Hyrax/Abif/Generate.hs b/src/Hyrax/Abif/Generate.hs --- a/src/Hyrax/Abif/Generate.hs +++ b/src/Hyrax/Abif/Generate.hs @@ -66,6 +66,20 @@ <> + +== Reverse reads + +A weighted FASTA can represent a reverse read. To do this add a `R` suffix to the weight. +The data you enter should be entered as if it was a forward read. This data will be complemented +and reversed before writing to the ABIF + +E.g. + +@ +> 1R +ACAG +@ + See README.md for additional details and examples -} module Hyrax.Abif.Generate @@ -74,6 +88,7 @@ , readWeightedFasta , iupac , unIupac + , complementNucleotides ) where import Protolude @@ -112,6 +127,7 @@ traverse_ (\(name, ab1) -> BS.writeFile (dest Txt.unpack name <> ".ab1") $ BSL.toStrict ab1) ab1s +{-! SECTION< gen_generateAb1_fn !-} -- | 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) = @@ -120,14 +136,17 @@ valsPerBase = trValsPerBase tr generatedFastaLen = (Txt.length $ trFasta tr) +{-! SECTION< gen_generateAb1_peak !-} -- 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..] +{-! SECTION> gen_generateAb1_peak !-} -- Sample name (from the FASTA name) sampleName = fst . Txt.breakOn "_" $ fName +{-! SECTION< gen_generateAb1_abif !-} -- Create the ABIF directories dirs = [ mkData 9 $ trData09G tr -- G , mkData 10 $ trData10A tr -- A @@ -149,33 +168,48 @@ , aRootDir = mkRoot , aDirs = dirs } +{-! SECTION> gen_generateAb1_abif !-} in -- Generate the data B.runPut (putAbif abif) +{-! SECTION> gen_generateAb1_fn !-} +{-! SECTION< gen_generateTraceData !-} +{-! SECTION< gen_generateTraceData_type !-} -- | Generate the traces for the AB1 from the parsed weighted FASTA generateTraceData :: [(Double, Text)] -> TraceData generateTraceData weighted = +{-! SECTION> gen_generateTraceData_type !-} let +{-! SECTION< gen_generateTraceData_weighted !-} weightedNucs' = (\(w, ns) -> (w,) . unIupac <$> Txt.unpack ns) <$> weighted weightedNucs = Lst.transpose weightedNucs' +{-! SECTION> gen_generateTraceData_weighted !-} - -- Values for a base that was present. This defines the shape of the chromatogram curve, and defines the number of values per base +{-! SECTION< gen_generateTraceData_curve !-} + -- 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 +{-! SECTION> gen_generateTraceData_curve !-} +{-! SECTION< gen_generateTraceData_traces !-} -- 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 +{-! SECTION> gen_generateTraceData_traces !-} +{-! SECTION< gen_generateTraceData_fasta !-} -- Create fasta sequence for the trace fastaSeq = concat <$> (snd <<$>> weightedNucs) fasta = Txt.pack $ iupac fastaSeq +{-! SECTION> gen_generateTraceData_fasta !-} in +{-! SECTION< gen_generateTraceData_ret !-} TraceData { trData09G = data09G , trData10A = data10A , trData11T = data11T @@ -183,8 +217,10 @@ , trFasta = fasta , trValsPerBase = valsPerBase } +{-! SECTION> gen_generateTraceData_ret !-} where +{-! SECTION< gen_generateTraceData_getWeightedTrace !-} getWeightedTrace :: [Int] -> Char -> [(Double, [Char])] -> [Int16] getWeightedTrace curve nuc ws = let @@ -194,10 +230,13 @@ wave = floor . (score *) . fromIntegral <$> curve in wave +{-! SECTION> gen_generateTraceData_getWeightedTrace !-} +{-! SECTION> gen_generateTraceData !-} --- | 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 +-- | Read a weighted FASTA file. See the module documentation for details on the format of the weighted FASTA +-- Reads with a weight followed by an `R` are reverse reads, and the AB1 generated will contain the complemeted +-- sequence. -- -- e.g. weighted FASTA -- @@ -214,14 +253,15 @@ -- The result data has the type -- -- @ --- ('Text', [('Double', 'Text')]) --- ^ ^ ^ --- | | | --- file name -------------+ | +---- read --- | --- +---- weight +-- [('Double', 'Text')] +-- ^ ^ +-- | | +-- | +---- read +-- | +-- +---- weight -- @ -- +{-! SECTION< gen_readWeightedFasta !-} readWeightedFasta :: ByteString -> Either Text [(Double, Text)] readWeightedFasta fastaData = case parseFasta $ TxtE.decodeUtf8 fastaData of @@ -236,14 +276,36 @@ Right r -> Right r readWeighted :: Fasta -> Either Text (Double, Text) - readWeighted (Fasta hdr dta) = + readWeighted (Fasta hdr' dta) = + let (processNucs, hdr) = + -- If there is a 'R' suffix, then generate a reverse sequence + -- Which means complement each nucleotide and then reverse the string + if Txt.isSuffixOf "R" hdr' + then (Txt.reverse . complementNucleotides, Txt.strip . Txt.dropEnd 1 $ hdr') + else (identity, hdr') + in + case (readMaybe . Txt.unpack $ hdr :: Maybe Double) of - Just weight -> Right (min 1 . max 0 $ weight, Txt.strip dta) + Just weight -> Right (min 1 . max 0 $ weight, processNucs $ Txt.strip dta) Nothing -> Left $ "Invalid header reading, expecting numeric weight, got: " <> hdr +{-! SECTION> gen_readWeightedFasta !-} +{-! SECTION< gen_readWeightedFastas !-} -- | Read all FASTA files in a directory +-- +-- The result data has the type +-- +-- @ +-- [ ('Text', [('Double', 'Text')]) ] +-- ^ ^ ^ +-- | | | +-- file name -------------+ | +---- read +-- | +-- +---- weight +-- @ +-- readWeightedFastas :: FilePath -> IO (Either Text [(Text, [(Double, Text)])]) readWeightedFastas source = do files <- filter (Txt.isSuffixOf ".fasta" . Txt.pack) <$> getFiles source @@ -253,6 +315,7 @@ case sequenceA $ readWeightedFasta <$> contents of Left e -> pure . Left $ e Right rs -> pure . Right $ zip names rs +{-! SECTION> gen_readWeightedFastas !-} -- | Find all files in a directory @@ -262,8 +325,11 @@ filterM Dir.doesFileExist entries +{-! SECTION< gen_unIupac_fn !-} -- | Convert a IUPAC ambiguity code to the set of nucleotides it represents +{-! SECTION< gen_unIupac_type !-} unIupac :: Char -> [Char] +{-! SECTION> gen_unIupac_type !-} unIupac c = case c of 'T' -> "T" @@ -286,6 +352,7 @@ 'X' -> "GATC" _ -> "" +{-! SECTION> gen_unIupac_fn !-} -- | Given a set of nucleotides get the IUPAC ambiguity code @@ -318,3 +385,23 @@ (False, True, True, True ) -> 'B' (True, True, True, True ) -> 'N' _ -> '_' + + +{-! SECTION< gen_complement !-} +-- | Return the complement of a nucelotide string +complementNucleotides :: Text -> Text +complementNucleotides ns = + let + un = unIupac <$> Txt.unpack ns + comp = complementNuc <<$>> un + iu = iupac comp + in + Txt.pack iu + + where + complementNuc 'A' = 'T' + complementNuc 'G' = 'C' + complementNuc 'T' = 'A' + complementNuc 'C' = 'G' + complementNuc x = x +{-! SECTION> gen_complement !-} diff --git a/src/Hyrax/Abif/Read.hs b/src/Hyrax/Abif/Read.hs --- a/src/Hyrax/Abif/Read.hs +++ b/src/Hyrax/Abif/Read.hs @@ -47,11 +47,14 @@ import Hyrax.Abif +{-! SECTION< read_readAbif !-} -- | Read and parse an AB1 file readAbif :: FilePath -> IO (Either Text Abif) readAbif path = getAbif <$> BSL.readFile path +{-! SECTION> read_readAbif !-} +{-! SECTION< read_getAbif !-} -- | Parse an AB1 from a 'ByteString' getAbif :: BSL.ByteString -> Either Text Abif getAbif bs = do @@ -63,23 +66,31 @@ 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) + Left (_, _, e) -> Left ("Error reading " + <> show (dElemNum rootDir) + <> " directories (at " <> show (dDataOffset rootDir) <> "): " + <> Txt.pack e + ) pure $ Abif header rootDir ds +{-! SECTION> read_getAbif !-} +{-! SECTION< read_clear !-} -- | Removes all data from the ABIF's directories clearAbif :: Abif -> Abif clearAbif a = a { aRootDir = clear $ aRootDir a - , aDirs = clear <$> aDirs 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 = "" } +{-! SECTION> read_clear !-} +{-! SECTION< read_getDebug !-} -- | 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 @@ -99,6 +110,7 @@ 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] } +{-! SECTION> read_getDebug !-} y -> -- For non-array entries @@ -175,6 +187,7 @@ pure (c:cs) +{-! SECTION< read_getStrings !-} -- | Parse a 'ElemPString' getPString :: B.Get Text getPString = do @@ -186,23 +199,29 @@ getCString :: Int -> B.Get Text getCString sz = TxtE.decodeUtf8 <$> B.getByteString (sz - 1) +{-! SECTION> read_getStrings !-} +{-! SECTION< read_getHeader !-} -- | Parse the ABIF 'Header' getHeader :: B.Get Header getHeader = Header <$> (TxtE.decodeUtf8 <$> B.getByteString 4) <*> (fromIntegral <$> B.getInt16be) +{-! SECTION> read_getHeader !-} +{-! SECTION< read_getRoot !-} -- | Parse the root ('Header' and 'Directory') getRoot :: BSL.ByteString -> B.Get (Header, Directory) getRoot bs = do h <- getHeader rd <- getDirectory bs pure (h, rd) +{-! SECTION> read_getRoot !-} +{-! SECTION< read_getDirectory !-} -- | Parse a single 'Directory' entry and read its data getDirectory :: BSL.ByteString -> B.Get Directory getDirectory bs = do @@ -221,7 +240,14 @@ 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 + 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 @@ -236,8 +262,10 @@ , dData = dataBytes , dDataDebug = [] } +{-! SECTION> read_getDirectory !-} +{-! SECTION< read_getDirectories !-} -- | Parse all the directoy entries getDirectories :: BSL.ByteString -> [Directory] -> Int -> B.Get [Directory] getDirectories _ acc 0 = pure acc @@ -246,4 +274,4 @@ B.skip 4 -- Skip the reserved field getDirectories bs (acc <> [d]) (more - 1) - +{-! SECTION> read_getDirectories !-} diff --git a/src/Hyrax/Abif/Write.hs b/src/Hyrax/Abif/Write.hs --- a/src/Hyrax/Abif/Write.hs +++ b/src/Hyrax/Abif/Write.hs @@ -44,10 +44,13 @@ import Hyrax.Abif +{-! SECTION< write_base !-} -- | Used to specify the base order for the FWO directry entry, see 'mkBaseOrder' data Base = BaseA | BaseC | BaseG | BaseT +{-! SECTION> write_base !-} +{-! SECTION< write_create !-} -- | Write an 'Abif' to a 'ByteString' createAbifBytes :: Abif -> BSL.ByteString createAbifBytes ab1 = @@ -59,19 +62,21 @@ writeAbif destPath ab1 = do let b = createAbifBytes ab1 BS.writeFile destPath $ BSL.toStrict b +{-! SECTION> write_create !-} +{-! SECTION< write_putAbif !-} -- | 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 + -- Data starts at offset 128 + let startDataOffset = 128 -- Write the root directory entry putDirectory (startDataOffset + dataSize) $ root { dDataSize = 28 * length dirs , dElemNum = length dirs @@ -91,8 +96,10 @@ pure $ if dDataSize dir > 4 then offset + dDataSize dir else offset +{-! SECTION> write_putAbif !-} +{-! SECTION< write_strings !-} -- | Write 'Text' putTextStr :: Text -> B.Put putTextStr t = B.putByteString $ TxtE.encodeUtf8 t @@ -103,15 +110,19 @@ putPStr t = do B.putInt8 . fromIntegral $ Txt.length t B.putByteString $ TxtE.encodeUtf8 t +{-! SECTION> write_strings !-} +{-! SECTION< write_putHeader !-} -- | Write a 'Header' putHeader :: Header -> B.Put putHeader h = do putTextStr $ hName h B.putInt16be . fromIntegral $ hVersion h +{-! SECTION> write_putHeader !-} +{-! SECTION< write_putDirectory !-} -- | Write a 'Directory' putDirectory :: Int -> Directory -> B.Put putDirectory dirOffset d = do @@ -129,6 +140,7 @@ else B.putLazyByteString . BSL.take 4 $ dData d <> "\0\0\0\0" B.putInt32be 0 -- reserved / datahandle +{-! SECTION> write_putDirectory !-} -- | Create a 'Header' @@ -192,6 +204,7 @@ , dDataSize = fromIntegral (BSL.length sampleName) } +{-! SECTION< write_mk !-} -- | Create a base order (FWO_) 'Directory' entry data mkBaseOrder :: Base -> Base -> Base -> Base -> Directory mkBaseOrder w x y z = @@ -229,6 +242,7 @@ , dData = B.runPut $ B.putInt16be lane , dDataDebug = [] } +{-! SECTION> write_mk !-} -- | Create a called bases (PBAS) 'Directory' entry and data @@ -328,7 +342,9 @@ , dElemNum = length ds } +{-! SECTION< write_addDirectory !-} -- | Add a directory to an 'Abif' addDirectory :: Abif -> Directory -> Abif addDirectory abif dir = abif { aDirs = aDirs abif <> [dir] } +{-! SECTION> write_addDirectory !-} diff --git a/test/AbifTests.hs b/test/AbifTests.hs --- a/test/AbifTests.hs +++ b/test/AbifTests.hs @@ -20,6 +20,19 @@ import qualified Hyrax.Abif.Generate as H import Generators + +-- | Test that complementing nucleotides is reversible +prop_complementNucs :: Property +prop_complementNucs = property $ do + nucs' <- forAll $ nucsGen + let nucs = Txt.replace "X" "N" nucs' + + let n1 = H.complementNucleotides nucs + let n2 = Txt.replace "X" "N" $ H.complementNucleotides n1 + + nucs === n2 + + -- | Test that an ab1 (write, read, write, read) results in the original data prop_roundtrip :: Property prop_roundtrip = property $ do @@ -112,6 +125,7 @@ c <- getFn cs <- readArray getFn pure (c:cs) + tests :: IO Bool