diff --git a/CHANGELOG.md b/CHANGELOG.md index 3d41825..7555a09 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,9 +6,16 @@ The format is based on Keep a Changelog and this project adheres to Semantic Ver ## [UNRELEASED](https://github.com/camilogarciabotero/GeneFinder.jl/compare/v0.0.10...main) -- Add `write_gff` methods. -- Clean up comments and docs. -- Add `iscoding` method. + +## [0.2.0] + +- A more stable ORF type that asserts for incorrect ORFs. +- Extended `Base.getindex` method to better handle ORF enabling: `seq[orf]` regardless of the ORF strand. +- The `write_*` API is much more stable and its been rewriten to handle IOStreams and `IOBuffers` as well. +- The new `record_orfs_fna` and `record_orfs_faa` methods intearact with the FASTX interface. +- Codebase has been reduced. Direct String manipulation is not enabled for `findorfs` and `get_*` methods. +- Docstrings have been updated following method updates. +- Better tests and now `Aqua.jl` is used for several tests as well. ## [0.1.0] diff --git a/Manifest.toml b/Manifest.toml index 50d74a0..239b166 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,10 +2,7 @@ julia_version = "1.10.0" manifest_format = "2.0" -project_hash = "5f5bec69be3c8dc8bff8fce86338a1ef1e8ac046" - -[[deps.Artifacts]] -uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" +project_hash = "4f6774eaeaec568c8cb7a865a84c0c7d9a180727" [[deps.Automa]] deps = ["PrecompileTools", "TranscodingStreams"] @@ -13,29 +10,12 @@ git-tree-sha1 = "0da671c730d79b8f9a88a391556ec695ea921040" uuid = "67c07d97-cdcb-5c2c-af73-a7f9c32a568b" version = "1.0.2" -[[deps.Base64]] -uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" - [[deps.BioGenerics]] deps = ["TranscodingStreams"] git-tree-sha1 = "7bbc085aebc6faa615740b63756e4986c9e85a70" uuid = "47718e42-2ac5-11e9-14af-e5595289c2ea" version = "0.1.4" -[[deps.BioMarkovChains]] -deps = ["BioSequences", "PrecompileTools", "VectorizedKmers"] -git-tree-sha1 = "44da2d71c2e9a68babb324c4df69d87b9ff7a835" -uuid = "f861b655-cb5f-42ce-b66a-341b542d4f2c" -version = "0.9.0" - - [deps.BioMarkovChains.extensions] - DiscreteMarkovChainsExt = "DiscreteMarkovChains" - MarkovChainHammerExt = "MarkovChainHammer" - - [deps.BioMarkovChains.weakdeps] - DiscreteMarkovChains = "8abcb7ef-b365-4f7b-ac38-56893fb62f9f" - MarkovChainHammer = "38c40fd0-bccb-4723-b30d-b2caea0ad8d9" - [[deps.BioSequences]] deps = ["BioSymbols", "PrecompileTools", "Random", "Twiddle"] git-tree-sha1 = "6fdba8b4279460fef5674e9aa2dac7ef5be361d5" @@ -48,42 +28,10 @@ git-tree-sha1 = "e32a61f028b823a172c75e26865637249bb30dff" uuid = "3c28c6f8-a34d-59c4-9654-267d177fcfa9" version = "5.1.3" -[[deps.Compat]] -deps = ["UUIDs"] -git-tree-sha1 = "886826d76ea9e72b35fcd000e535588f7b60f21d" -uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" -version = "4.10.1" -weakdeps = ["Dates", "LinearAlgebra"] - - [deps.Compat.extensions] - CompatLinearAlgebraExt = "LinearAlgebra" - -[[deps.CompilerSupportLibraries_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" -version = "1.0.5+1" - -[[deps.DataAPI]] -git-tree-sha1 = "8da84edb865b0b5b0100c0666a9bc9a0b71c553c" -uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a" -version = "1.15.0" - -[[deps.DataStructures]] -deps = ["Compat", "InteractiveUtils", "OrderedCollections"] -git-tree-sha1 = "3dbd312d370723b6bb43ba9d02fc36abade4518d" -uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" -version = "0.18.15" - [[deps.Dates]] deps = ["Printf"] uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" -[[deps.DocStringExtensions]] -deps = ["LibGit2"] -git-tree-sha1 = "2fb1e02f2b635d0845df5d7c167fec4dd739b00d" -uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -version = "0.9.3" - [[deps.FASTX]] deps = ["Automa", "BioGenerics", "PrecompileTools", "StringViews", "TranscodingStreams"] git-tree-sha1 = "bff5d62bf5e1c382a370ac701bcaea9a24115ac6" @@ -94,85 +42,10 @@ weakdeps = ["BioSequences"] [deps.FASTX.extensions] BioSequencesExt = "BioSequences" -[[deps.InteractiveUtils]] -deps = ["Markdown"] -uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" - -[[deps.IrrationalConstants]] -git-tree-sha1 = "630b497eafcc20001bba38a4651b327dcfc491d2" -uuid = "92d709cd-6900-40b7-9082-c6be49f344b6" -version = "0.2.2" - [[deps.IterTools]] -git-tree-sha1 = "274c38bd733f9d29036d0a73658fff1dc1d3a065" +git-tree-sha1 = "42d5f897009e7ff2cf88db414a389e5ed1bdd023" uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e" -version = "1.9.0" - -[[deps.LibGit2]] -deps = ["Base64", "LibGit2_jll", "NetworkOptions", "Printf", "SHA"] -uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" - -[[deps.LibGit2_jll]] -deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll"] -uuid = "e37daf67-58a4-590a-8e99-b0245dd2ffc5" -version = "1.6.4+0" - -[[deps.LibSSH2_jll]] -deps = ["Artifacts", "Libdl", "MbedTLS_jll"] -uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8" -version = "1.11.0+1" - -[[deps.Libdl]] -uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" - -[[deps.LinearAlgebra]] -deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"] -uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" - -[[deps.LogExpFunctions]] -deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"] -git-tree-sha1 = "7d6dd4e9212aebaeed356de34ccf262a3cd415aa" -uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688" -version = "0.3.26" - - [deps.LogExpFunctions.extensions] - LogExpFunctionsChainRulesCoreExt = "ChainRulesCore" - LogExpFunctionsChangesOfVariablesExt = "ChangesOfVariables" - LogExpFunctionsInverseFunctionsExt = "InverseFunctions" - - [deps.LogExpFunctions.weakdeps] - ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" - ChangesOfVariables = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0" - InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112" - -[[deps.Markdown]] -deps = ["Base64"] -uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" - -[[deps.MbedTLS_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" -version = "2.28.2+1" - -[[deps.Missings]] -deps = ["DataAPI"] -git-tree-sha1 = "f66bdc5de519e8f8ae43bdc598782d35a25b1272" -uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28" -version = "1.1.0" - -[[deps.NetworkOptions]] -uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" -version = "1.2.0" - -[[deps.OpenBLAS_jll]] -deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] -uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" -version = "0.3.23+2" - -[[deps.OrderedCollections]] -git-tree-sha1 = "dfdf5519f235516220579f949664f1bf44e741c5" -uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" -version = "1.6.3" +version = "1.10.0" [[deps.PrecompileTools]] deps = ["Preferences"] @@ -198,47 +71,11 @@ uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" version = "0.7.0" -[[deps.Serialization]] -uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" - -[[deps.SortingAlgorithms]] -deps = ["DataStructures"] -git-tree-sha1 = "66e0a8e672a0bdfca2c3f5937efb8538b9ddc085" -uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c" -version = "1.2.1" - -[[deps.SparseArrays]] -deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"] -uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" -version = "1.10.0" - -[[deps.Statistics]] -deps = ["LinearAlgebra", "SparseArrays"] -uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" -version = "1.10.0" - -[[deps.StatsAPI]] -deps = ["LinearAlgebra"] -git-tree-sha1 = "1ff449ad350c9c4cbc756624d6f8a8c3ef56d3ed" -uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0" -version = "1.7.0" - -[[deps.StatsBase]] -deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"] -git-tree-sha1 = "1d77abd07f617c4868c33d4f5b9e1dbb2643c9cf" -uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" -version = "0.34.2" - [[deps.StringViews]] git-tree-sha1 = "f7b06677eae2571c888fd686ba88047d8738b0e3" uuid = "354b36f9-a18e-4713-926e-db85100087ba" version = "1.3.3" -[[deps.SuiteSparse_jll]] -deps = ["Artifacts", "Libdl", "libblastrampoline_jll"] -uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c" -version = "7.2.1+1" - [[deps.TOML]] deps = ["Dates"] uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" @@ -261,23 +98,5 @@ git-tree-sha1 = "29509c4862bfb5da9e76eb6937125ab93986270a" uuid = "7200193e-83a8-5a55-b20d-5d36d44a0795" version = "1.1.2" -[[deps.UUIDs]] -deps = ["Random", "SHA"] -uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" - [[deps.Unicode]] uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" - -[[deps.VectorizedKmers]] -git-tree-sha1 = "eaf497f1b720ae17dc52463dfe5fed5f74b0e4df" -uuid = "2ef45bd6-4c2b-43d8-88b3-40597287359a" -version = "0.8.1" -weakdeps = ["BioSequences"] - - [deps.VectorizedKmers.extensions] - BioSequencesExt = "BioSequences" - -[[deps.libblastrampoline_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" -version = "5.8.0+1" diff --git a/Project.toml b/Project.toml index 1209169..093d413 100644 --- a/Project.toml +++ b/Project.toml @@ -1,27 +1,24 @@ name = "GeneFinder" uuid = "2bc6ee46-2158-4f5a-a720-91cb2d7cee64" authors = ["Camilo García-Botero"] -version = "0.1.0" +version = "0.2.0" [deps] -BioMarkovChains = "f861b655-cb5f-42ce-b66a-341b542d4f2c" BioSequences = "7e6ae17a-c86d-528c-b3b9-7f778a29fe59" FASTX = "c2308a5c-f048-11e8-3e8a-31650f418d12" IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" -StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" [compat] -BioMarkovChains = "0.9.0" BioSequences = "3" FASTX = "2" IterTools = "1.4" PrecompileTools = "1" -StatsBase = "0.34.0" julia = "1" [extras] +Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test"] +test = ["Test", "Aqua"] diff --git a/README.md b/README.md index 86b6899..fe1a30a 100644 --- a/README.md +++ b/README.md @@ -12,6 +12,7 @@ [![License](https://img.shields.io/badge/license-MIT-green.svg)](https://github.com/camilogarciabotero/GeneFinder.jl/blob/main/LICENSE) [![Repo Status](https://www.repostatus.org/badges/latest/wip.svg)](https://www.repostatus.org/#wip) [![Downloads](https://shields.io/endpoint?url=https://pkgs.genieframework.com/api/v1/badge/GeneFinder&label=downloads)](https://pkgs.genieframework.com?packages=GeneFinder) +[![Aqua QA](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl) @@ -36,7 +37,7 @@ add GeneFinder If you are interested in the cutting edge of the development, please check out the master branch to try new features before release. -## Finding complete and internal (overlapped) ORFs +## Finding complete and overlapped ORFs The first implemented function is `findorfs` a very non-restrictive ORF finder function that will catch all ORFs in a dedicated structure. Note that this will catch random ORFs not necesarily genes since it has no ORFs size or overlapping condition contraints. Thus it might consider `aa"M*"` a posible encoding protein from the resulting ORFs. @@ -66,10 +67,10 @@ findorfs(seq) ORF(695:706, '+', 2) ``` -Two other functions (`getorfdna` and `getorfaa`) pass the sequence to `findorfs` take the ORFs and act as generators of the sequence, so this way the can be `collect`ed in the REPL as an standard output or writteen into a file more conviniently using the `FASTX` IO system: +Two other functions (`get_orfs_dna` and `get_orfs_aa`) are implemented to get the ORFs in DNA and amino acid sequences, respectively. They use the `findorfs` function to first get the ORFs and then get the correspondance array of `BioSequence` objects. ```julia -getorfdna(seq) +get_orfs_dna(seq) 12-element Vector{LongSubSeq{DNAAlphabet{4}}}: ATGCAACCCTGA @@ -86,78 +87,68 @@ getorfdna(seq) ATGCAACCCTGA ``` +## Writting ORF information into bioinformatic formats + +This package facilitates now the creation of `FASTA`, `BED`, and `GFF` files, specifically extracting Open Reading Frame (ORF) information from `BioSequence` instances, particularly those of type `NucleicSeqOrView{A} where A`, and then writing the information into the desired format. + +Functionality: + +The package provides four distinct functions for writing files in different formats: + +| Function | Description | +|-------------------|--------------------------------------------------------| +| `write_orfs_fna` | Writes nucleotide sequences in FASTA format. | +| `write_orfs_faa` | Writes amino acid sequences in FASTA format. | +| `write_orfs_bed` | Outputs information in BED format. | +| `write_orfs_gff` | Generates files in GFF format. | + + +All these functions support processing both `BioSequence` instances and external `FASTA` files. In the case of a `BioSequence` instace into external files, simply provide the path to the `FASTA` file using a `String` to the path. To demonstrate the use of the `write_*` methods with a `BioSequence`, consider the following example: + ```julia -getorfaa(seq) - -12-element Vector{LongSubSeq{AminoAcidAlphabet}}: - MQP* - MR* - MRRMAR* - MAR* - M* - MCPTAV* - MQP* - MHWLVLSI* - MSPHKAM* - M* - MCPTAA* - MQP* +using BioSequences, GeneFinder + +# > 180195.SAMN03785337.LFLS01000089 -> finds only 1 gene in Prodigal (from Pyrodigal tests) +seq = dna"AACCAGGGCAATATCAGTACCGCGGGCAATGCAACCCTGACTGCCGGCGGTAACCTGAACAGCACTGGCAATCTGACTGTGGGCGGTGTTACCAACGGCACTGCTACTACTGGCAACATCGCACTGACCGGTAACAATGCGCTGAGCGGTCCGGTCAATCTGAATGCGTCGAATGGCACGGTGACCTTGAACACGACCGGCAATACCACGCTCGGTAACGTGACGGCACAAGGCAATGTGACGACCAATGTGTCCAACGGCAGTCTGACGGTTACCGGCAATACGACAGGTGCCAACACCAACCTCAGTGCCAGCGGCAACCTGACCGTGGGTAACCAGGGCAATATCAGTACCGCAGGCAATGCAACCCTGACGGCCGGCGACAACCTGACGAGCACTGGCAATCTGACTGTGGGCGGCGTCACCAACGGCACGGCCACCACCGGCAACATCGCGCTGACCGGTAACAATGCACTGGCTGGTCCTGTCAATCTGAACGCGCCGAACGGCACCGTGACCCTGAACACAACCGGCAATACCACGCTGGGTAATGTCACCGCACAAGGCAATGTGACGACTAATGTGTCCAACGGCAGCCTGACAGTCGCTGGCAATACCACAGGTGCCAACACCAACCTGAGTGCCAGCGGCAATCTGACCGTGGGCAACCAGGGCAATATCAGTACCGCGGGCAATGCAACCCTGACTGCCGGCGGTAACCTGAGC" ``` -### Writting cds, proteins fastas, bed and gffs whether from a `LongSeq` or from a external fasta file. +Once a `BioSequence` object has been instantiated, the `write_orfs_fna` function proves useful for generating a `FASTA` file containing the nucleotide sequences of the ORFs. Notably, the `write_orfs*` methods support either an `IOStream` or an `IOBuffer` as an output argument, allowing flexibility in directing the output either to a file or a buffer. In the following example, we demonstrate writing the output directly to a file. ```julia -write_cds("cds.fasta", seq) +outfile = "LFLS01000089.fna" + +open(outfile, "w") do io + write_orfs_fna(seq, io) +end ``` ```bash -cat cds.fasta +cat LFLS01000089.fna ->location=29:40 strand=+ frame=2 +>ORF01 id=01 start=29 stop=40 strand=+ frame=2 ATGCAACCCTGA ->location=137:145 strand=+ frame=2 +>ORF02 id=02 start=137 stop=145 strand=+ frame=2 ATGCGCTGA ->location=164:184 strand=+ frame=2 +>ORF03 id=03 start=164 stop=184 strand=+ frame=2 ATGCGTCGAATGGCACGGTGA ->location=173:184 strand=+ frame=2 +>ORF04 id=04 start=173 stop=184 strand=+ frame=2 ATGGCACGGTGA ->location=236:241 strand=+ frame=2 +>ORF05 id=05 start=236 stop=241 strand=+ frame=2 ATGTGA ->location=248:268 strand=+ frame=2 +>ORF06 id=06 start=248 stop=268 strand=+ frame=2 ATGTGTCCAACGGCAGTCTGA ->location=362:373 strand=+ frame=2 +>ORF07 id=07 start=362 stop=373 strand=+ frame=2 ATGCAACCCTGA ->location=470:496 strand=+ frame=2 +>ORF08 id=08 start=470 stop=496 strand=+ frame=2 ATGCACTGGCTGGTCCTGTCAATCTGA ->location=551:574 strand=+ frame=2 +>ORF09 id=09 start=551 stop=574 strand=+ frame=2 ATGTCACCGCACAAGGCAATGTGA ->location=569:574 strand=+ frame=2 +>ORF10 id=10 start=569 stop=574 strand=+ frame=2 ATGTGA ->location=581:601 strand=+ frame=2 +>ORF11 id=11 start=581 stop=601 strand=+ frame=2 ATGTGTCCAACGGCAGCCTGA ->location=695:706 strand=+ frame=2 +>ORF12 id=12 start=695 stop=706 strand=+ frame=2 ATGCAACCCTGA ``` -### Combining `FASTX` for reading and writing fastas - -```julia -using FASTX - -write_proteins("test/data/NC_001884.fasta", "proteins.fasta") -``` - -```bash -head proteins.fasta - ->location=41:145 strand=- frame=2 -MTQKRKGPIPAQFEITPILRFNFIFDLTATNSFH* ->location=41:172 strand=- frame=2 -MVLKDVIVNMTQKRKGPIPAQFEITPILRFNFIFDLTATNSFH* ->location=41:454 strand=- frame=2 -MSEHLSQKEKELKNKENFIFDKYESGIYSDELFLKRKAALDEEFKELQNAKNELNGLQDTQSEIDSNTVRNNINKIIDQYHIESSSEKKNELLRMVLKDVIVNMTQKRKGPIPAQFEITPILRFNFIFDLTATNSFH* ->location=41:472 strand=- frame=2 -MKTKKQMSEHLSQKEKELKNKENFIFDKYESGIYSDELFLKRKAALDEEFKELQNAKNELNGLQDTQSEIDSNTVRNNINKIIDQYHIESSSEKKNELLRMVLKDVIVNMTQKRKGPIPAQFEITPILRFNFIFDLTATNSFH* ->location=41:505 strand=- frame=2 -MLSKYEDDNSNMKTKKQMSEHLSQKEKELKNKENFIFDKYESGIYSDELFLKRKAALDEEFKELQNAKNELNGLQDTQSEIDSNTVRNNINKIIDQYHIESSSEKKNELLRMVLKDVIVNMTQKRKGPIPAQFEITPILRFNFIFDLTATNSFH* -``` \ No newline at end of file +This could also be done to writting a `FASTA` file with the nucleotide sequences of the ORFs using the `write_orfs_fna` function. Similarly for the `BED` and `GFF` files using the `write_orfs_bed` and `write_orfs_gff` functions respectively. \ No newline at end of file diff --git a/docs/Manifest.toml b/docs/Manifest.toml index e8874f3..d99e7c6 100644 --- a/docs/Manifest.toml +++ b/docs/Manifest.toml @@ -1,8 +1,8 @@ # This file is machine-generated - editing it directly is not advised -julia_version = "1.9.3" +julia_version = "1.10.0" manifest_format = "2.0" -project_hash = "6d443a9953f0c7dc7584294d2ada98dae33c8c16" +project_hash = "e796693959137d83427b789f5054fa800312cae0" [[deps.ANSIColoredPrinters]] git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c" @@ -23,24 +23,38 @@ uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" [[deps.Automa]] deps = ["PrecompileTools", "TranscodingStreams"] -git-tree-sha1 = "9dca1d1632c4f7fb320b8a766feb6b8290695f7c" +git-tree-sha1 = "0da671c730d79b8f9a88a391556ec695ea921040" uuid = "67c07d97-cdcb-5c2c-af73-a7f9c32a568b" -version = "1.0.1" +version = "1.0.2" [[deps.Base64]] uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" [[deps.BioGenerics]] deps = ["TranscodingStreams"] -git-tree-sha1 = "0b581906418b93231d391b5dd78831fdc2da0c82" +git-tree-sha1 = "7bbc085aebc6faa615740b63756e4986c9e85a70" uuid = "47718e42-2ac5-11e9-14af-e5595289c2ea" -version = "0.1.2" +version = "0.1.4" + +[[deps.BioMarkovChains]] +deps = ["BioSequences", "PrecompileTools", "VectorizedKmers"] +git-tree-sha1 = "44da2d71c2e9a68babb324c4df69d87b9ff7a835" +uuid = "f861b655-cb5f-42ce-b66a-341b542d4f2c" +version = "0.9.0" + + [deps.BioMarkovChains.extensions] + DiscreteMarkovChainsExt = "DiscreteMarkovChains" + MarkovChainHammerExt = "MarkovChainHammer" + + [deps.BioMarkovChains.weakdeps] + DiscreteMarkovChains = "8abcb7ef-b365-4f7b-ac38-56893fb62f9f" + MarkovChainHammer = "38c40fd0-bccb-4723-b30d-b2caea0ad8d9" [[deps.BioSequences]] deps = ["BioSymbols", "PrecompileTools", "Random", "Twiddle"] -git-tree-sha1 = "6cf406ea0d5bc901eabe55fca609007bce2db026" +git-tree-sha1 = "6fdba8b4279460fef5674e9aa2dac7ef5be361d5" uuid = "7e6ae17a-c86d-528c-b3b9-7f778a29fe59" -version = "3.1.5" +version = "3.1.6" [[deps.BioSymbols]] deps = ["PrecompileTools"] @@ -49,47 +63,21 @@ uuid = "3c28c6f8-a34d-59c4-9654-267d177fcfa9" version = "5.1.3" [[deps.BitFlags]] -git-tree-sha1 = "43b1a4a8f797c1cddadf60499a8a077d4af2cd2d" +git-tree-sha1 = "2dc09997850d68179b69dafb58ae806167a32b1b" uuid = "d1d4a3ce-64b1-5f1a-9ba4-7e7e69966f35" -version = "0.1.7" +version = "0.1.8" [[deps.CodecZlib]] deps = ["TranscodingStreams", "Zlib_jll"] -git-tree-sha1 = "02aa26a4cf76381be7f66e020a3eddeb27b0a092" +git-tree-sha1 = "cd67fc487743b2f0fd4380d4cbd3a24660d0eec8" uuid = "944b1d66-785c-5afd-91f1-9de20f533193" -version = "0.7.2" - -[[deps.Compat]] -deps = ["UUIDs"] -git-tree-sha1 = "e460f044ca8b99be31d35fe54fc33a5c33dd8ed7" -uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" -version = "4.9.0" -weakdeps = ["Dates", "LinearAlgebra"] - - [deps.Compat.extensions] - CompatLinearAlgebraExt = "LinearAlgebra" - -[[deps.CompilerSupportLibraries_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" -version = "1.0.5+0" +version = "0.7.3" [[deps.ConcurrentUtilities]] deps = ["Serialization", "Sockets"] -git-tree-sha1 = "5372dbbf8f0bdb8c700db5367132925c0771ef7e" +git-tree-sha1 = "8cfa272e8bdedfa88b6aefbbca7c19f1befac519" uuid = "f0e56b4a-5159-44fe-b623-3e5288b988bb" -version = "2.2.1" - -[[deps.DataAPI]] -git-tree-sha1 = "8da84edb865b0b5b0100c0666a9bc9a0b71c553c" -uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a" -version = "1.15.0" - -[[deps.DataStructures]] -deps = ["Compat", "InteractiveUtils", "OrderedCollections"] -git-tree-sha1 = "3dbd312d370723b6bb43ba9d02fc36abade4518d" -uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" -version = "0.18.15" +version = "2.3.0" [[deps.Dates]] deps = ["Printf"] @@ -102,10 +90,10 @@ uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" version = "0.9.3" [[deps.Documenter]] -deps = ["ANSIColoredPrinters", "AbstractTrees", "Base64", "Dates", "DocStringExtensions", "Downloads", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "MarkdownAST", "Pkg", "PrecompileTools", "REPL", "RegistryInstances", "SHA", "Test", "Unicode"] -git-tree-sha1 = "43a76dfac223a3a0d7d33443f7c9154fe75bb264" +deps = ["ANSIColoredPrinters", "AbstractTrees", "Base64", "Dates", "DocStringExtensions", "Downloads", "Git", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "MarkdownAST", "Pkg", "PrecompileTools", "REPL", "RegistryInstances", "SHA", "Test", "Unicode"] +git-tree-sha1 = "2613dbec8f4748273bbe30ba71fd5cb369966bac" uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -version = "1.0.0" +version = "1.2.1" [[deps.Downloads]] deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"] @@ -114,15 +102,21 @@ version = "1.6.0" [[deps.ExceptionUnwrapping]] deps = ["Test"] -git-tree-sha1 = "e90caa41f5a86296e014e148ee061bd6c3edec96" +git-tree-sha1 = "dcb08a0d93ec0b1cdc4af184b26b591e9695423a" uuid = "460bff9d-24e4-43bc-9d9f-a8973cb893f4" -version = "0.1.9" +version = "0.1.10" + +[[deps.Expat_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "4558ab818dcceaab612d1bb8c19cee87eda2b83c" +uuid = "2e619515-83b5-522b-bb60-26c02a35a201" +version = "2.5.0+0" [[deps.FASTX]] deps = ["Automa", "BioGenerics", "PrecompileTools", "StringViews", "TranscodingStreams"] -git-tree-sha1 = "c05aede7ef67db392d9c779257e709de72504447" +git-tree-sha1 = "bff5d62bf5e1c382a370ac701bcaea9a24115ac6" uuid = "c2308a5c-f048-11e8-3e8a-31650f418d12" -version = "2.1.3" +version = "2.1.4" weakdeps = ["BioSequences"] [deps.FASTX.extensions] @@ -132,16 +126,28 @@ weakdeps = ["BioSequences"] uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee" [[deps.GeneFinder]] -deps = ["BioSequences", "FASTX", "IterTools", "PrecompileTools", "StatsBase", "TestItemRunner", "TestItems"] +deps = ["BioSequences", "FASTX", "IterTools", "PrecompileTools"] path = ".." uuid = "2bc6ee46-2158-4f5a-a720-91cb2d7cee64" -version = "0.0.23" +version = "0.1.0" + +[[deps.Git]] +deps = ["Git_jll"] +git-tree-sha1 = "51764e6c2e84c37055e846c516e9015b4a291c7d" +uuid = "d7ba0133-e1db-5d97-8f8c-041e4b3a1eb2" +version = "1.3.0" + +[[deps.Git_jll]] +deps = ["Artifacts", "Expat_jll", "JLLWrappers", "LibCURL_jll", "Libdl", "Libiconv_jll", "OpenSSL_jll", "PCRE2_jll", "Zlib_jll"] +git-tree-sha1 = "bb8f7cc77ec1152414b2af6db533d9471cfbb2d1" +uuid = "f8c6e375-362e-5223-8a59-34ff63f689eb" +version = "2.42.0+0" [[deps.HTTP]] deps = ["Base64", "CodecZlib", "ConcurrentUtilities", "Dates", "ExceptionUnwrapping", "Logging", "LoggingExtras", "MbedTLS", "NetworkOptions", "OpenSSL", "Random", "SimpleBufferStream", "Sockets", "URIs", "UUIDs"] -git-tree-sha1 = "19e974eced1768fb46fd6020171f2cec06b1edb5" +git-tree-sha1 = "abbbb9ec3afd783a7cbd82ef01dcd088ea051398" uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3" -version = "1.9.15" +version = "1.10.1" [[deps.IOCapture]] deps = ["Logging", "Random"] @@ -153,15 +159,10 @@ version = "0.2.3" deps = ["Markdown"] uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" -[[deps.IrrationalConstants]] -git-tree-sha1 = "630b497eafcc20001bba38a4651b327dcfc491d2" -uuid = "92d709cd-6900-40b7-9082-c6be49f344b6" -version = "0.2.2" - [[deps.IterTools]] -git-tree-sha1 = "4ced6667f9974fc5c5943fa5e2ef1ca43ea9e450" +git-tree-sha1 = "42d5f897009e7ff2cf88db414a389e5ed1bdd023" uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e" -version = "1.8.0" +version = "1.10.0" [[deps.JLLWrappers]] deps = ["Artifacts", "Preferences"] @@ -176,35 +177,42 @@ uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" version = "0.21.4" [[deps.LazilyInitializedFields]] -git-tree-sha1 = "410fe4739a4b092f2ffe36fcb0dcc3ab12648ce1" +git-tree-sha1 = "8f7f3cabab0fd1800699663533b6d5cb3fc0e612" uuid = "0e77f7df-68c5-4e49-93ce-4cd80f5598bf" -version = "1.2.1" +version = "1.2.2" [[deps.LibCURL]] deps = ["LibCURL_jll", "MozillaCACerts_jll"] uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21" -version = "0.6.3" +version = "0.6.4" [[deps.LibCURL_jll]] deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"] uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0" -version = "7.84.0+0" +version = "8.4.0+0" [[deps.LibGit2]] -deps = ["Base64", "NetworkOptions", "Printf", "SHA"] +deps = ["Base64", "LibGit2_jll", "NetworkOptions", "Printf", "SHA"] uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" +[[deps.LibGit2_jll]] +deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll"] +uuid = "e37daf67-58a4-590a-8e99-b0245dd2ffc5" +version = "1.6.4+0" + [[deps.LibSSH2_jll]] deps = ["Artifacts", "Libdl", "MbedTLS_jll"] uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8" -version = "1.10.2+0" +version = "1.11.0+1" [[deps.Libdl]] uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" -[[deps.LinearAlgebra]] -deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"] -uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +[[deps.Libiconv_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "f9557a255370125b405568f9767d6d195822a175" +uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531" +version = "1.17.0+0" [[deps.LiveServer]] deps = ["HTTP", "LoggingExtras", "MIMEs", "Pkg", "Sockets", "Test"] @@ -212,30 +220,14 @@ git-tree-sha1 = "24d05efe53436b22a42bf2ae459f47c48b0c2603" uuid = "16fef848-5104-11e9-1b77-fb7a48bbb589" version = "1.2.7" -[[deps.LogExpFunctions]] -deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"] -git-tree-sha1 = "7d6dd4e9212aebaeed356de34ccf262a3cd415aa" -uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688" -version = "0.3.26" - - [deps.LogExpFunctions.extensions] - LogExpFunctionsChainRulesCoreExt = "ChainRulesCore" - LogExpFunctionsChangesOfVariablesExt = "ChangesOfVariables" - LogExpFunctionsInverseFunctionsExt = "InverseFunctions" - - [deps.LogExpFunctions.weakdeps] - ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" - ChangesOfVariables = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0" - InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112" - [[deps.Logging]] uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" [[deps.LoggingExtras]] deps = ["Dates", "Logging"] -git-tree-sha1 = "0d097476b6c381ab7906460ef1ef1638fbce1d91" +git-tree-sha1 = "c1dd6d7978c12545b4179fb6153b9250c96b0075" uuid = "e6f89c97-d47a-5376-807f-9c37f3926c36" -version = "1.0.2" +version = "1.0.3" [[deps.MIMEs]] git-tree-sha1 = "65f28ad4b594aebe22157d6fac869786a255b7eb" @@ -248,43 +240,32 @@ uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" [[deps.MarkdownAST]] deps = ["AbstractTrees", "Markdown"] -git-tree-sha1 = "e8513266815200c0c8f522d6d44ffb5e9b366ae4" +git-tree-sha1 = "465a70f0fc7d443a00dcdc3267a497397b8a3899" uuid = "d0879d2d-cac2-40c8-9cee-1863dc0c7391" -version = "0.1.1" +version = "0.1.2" [[deps.MbedTLS]] -deps = ["Dates", "MbedTLS_jll", "MozillaCACerts_jll", "Random", "Sockets"] -git-tree-sha1 = "03a9b9718f5682ecb107ac9f7308991db4ce395b" +deps = ["Dates", "MbedTLS_jll", "MozillaCACerts_jll", "NetworkOptions", "Random", "Sockets"] +git-tree-sha1 = "c067a280ddc25f196b5e7df3877c6b226d390aaf" uuid = "739be429-bea8-5141-9913-cc70e7f3736d" -version = "1.1.7" +version = "1.1.9" [[deps.MbedTLS_jll]] deps = ["Artifacts", "Libdl"] uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" -version = "2.28.2+0" - -[[deps.Missings]] -deps = ["DataAPI"] -git-tree-sha1 = "f66bdc5de519e8f8ae43bdc598782d35a25b1272" -uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28" -version = "1.1.0" +version = "2.28.2+1" [[deps.Mmap]] uuid = "a63ad114-7e13-5084-954f-fe012c677804" [[deps.MozillaCACerts_jll]] uuid = "14a3606d-f60d-562e-9121-12d972cd8159" -version = "2022.10.11" +version = "2023.1.10" [[deps.NetworkOptions]] uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" version = "1.2.0" -[[deps.OpenBLAS_jll]] -deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] -uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" -version = "0.3.21+4" - [[deps.OpenSSL]] deps = ["BitFlags", "Dates", "MozillaCACerts_jll", "OpenSSL_jll", "Sockets"] git-tree-sha1 = "51901a49222b09e3743c65b8847687ae5fc78eb2" @@ -293,25 +274,25 @@ version = "1.4.1" [[deps.OpenSSL_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "e78db7bd5c26fc5a6911b50a47ee302219157ea8" +git-tree-sha1 = "cc6e1927ac521b659af340e0ca45828a3ffc748f" uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95" -version = "3.0.10+0" +version = "3.0.12+0" -[[deps.OrderedCollections]] -git-tree-sha1 = "2e73fe17cac3c62ad1aebe70d44c963c3cfdc3e3" -uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" -version = "1.6.2" +[[deps.PCRE2_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "efcefdf7-47ab-520b-bdef-62a2eaa19f15" +version = "10.42.0+1" [[deps.Parsers]] deps = ["Dates", "PrecompileTools", "UUIDs"] -git-tree-sha1 = "716e24b21538abc91f6205fd1d8363f39b442851" +git-tree-sha1 = "8489905bcdbcfac64d1daa51ca07c0d8f0283821" uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" -version = "2.7.2" +version = "2.8.1" [[deps.Pkg]] deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"] uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" -version = "1.9.2" +version = "1.10.0" [[deps.PrecompileTools]] deps = ["Preferences"] @@ -321,9 +302,9 @@ version = "1.2.0" [[deps.Preferences]] deps = ["TOML"] -git-tree-sha1 = "7eb1686b4f04b82f96ed7a4ea5890a4f0c7a09f1" +git-tree-sha1 = "00805cd429dcb4870060ff49ef443486c262e38e" uuid = "21216c6a-2e73-6563-6e65-726566657250" -version = "1.4.0" +version = "1.4.1" [[deps.Printf]] deps = ["Unicode"] @@ -334,7 +315,7 @@ deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"] uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" [[deps.Random]] -deps = ["SHA", "Serialization"] +deps = ["SHA"] uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" [[deps.RegistryInstances]] @@ -358,43 +339,11 @@ version = "1.1.0" [[deps.Sockets]] uuid = "6462fe0b-24de-5631-8697-dd941f90decc" -[[deps.SortingAlgorithms]] -deps = ["DataStructures"] -git-tree-sha1 = "c60ec5c62180f27efea3ba2908480f8055e17cee" -uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c" -version = "1.1.1" - -[[deps.SparseArrays]] -deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"] -uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" - -[[deps.Statistics]] -deps = ["LinearAlgebra", "SparseArrays"] -uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" -version = "1.9.0" - -[[deps.StatsAPI]] -deps = ["LinearAlgebra"] -git-tree-sha1 = "1ff449ad350c9c4cbc756624d6f8a8c3ef56d3ed" -uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0" -version = "1.7.0" - -[[deps.StatsBase]] -deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"] -git-tree-sha1 = "75ebe04c5bed70b91614d684259b661c9e6274a4" -uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" -version = "0.34.0" - [[deps.StringViews]] git-tree-sha1 = "f7b06677eae2571c888fd686ba88047d8738b0e3" uuid = "354b36f9-a18e-4713-926e-db85100087ba" version = "1.3.3" -[[deps.SuiteSparse_jll]] -deps = ["Artifacts", "Libdl", "Pkg", "libblastrampoline_jll"] -uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c" -version = "5.10.1+6" - [[deps.TOML]] deps = ["Dates"] uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" @@ -409,22 +358,14 @@ version = "1.10.0" deps = ["InteractiveUtils", "Logging", "Random", "Serialization"] uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -[[deps.TestItemRunner]] -deps = ["Pkg", "TOML", "Test", "TestItems", "UUIDs"] -git-tree-sha1 = "7351f1f7389e048ce974b2c176513b817c94b15f" -uuid = "f8b46487-2199-4994-9208-9a1283c18c0a" -version = "0.2.2" - -[[deps.TestItems]] -git-tree-sha1 = "8621ba2637b49748e2dc43ba3d84340be2938022" -uuid = "1c621080-faea-4a02-84b6-bbd5e436b8fe" -version = "0.1.1" - [[deps.TranscodingStreams]] -deps = ["Random", "Test"] -git-tree-sha1 = "9a6ae7ed916312b41236fcef7e0af564ef934769" +git-tree-sha1 = "1fbeaaca45801b4ba17c251dd8603ef24801dd84" uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa" -version = "0.9.13" +version = "0.10.2" +weakdeps = ["Random", "Test"] + + [deps.TranscodingStreams.extensions] + TestExt = ["Test", "Random"] [[deps.Twiddle]] git-tree-sha1 = "29509c4862bfb5da9e76eb6937125ab93986270a" @@ -432,9 +373,9 @@ uuid = "7200193e-83a8-5a55-b20d-5d36d44a0795" version = "1.1.2" [[deps.URIs]] -git-tree-sha1 = "b7a5e99f24892b6824a954199a45e9ffcc1c70f0" +git-tree-sha1 = "67db6cc7b3821e19ebe75791a9dd19c9b1188f2b" uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4" -version = "1.5.0" +version = "1.5.1" [[deps.UUIDs]] deps = ["Random", "SHA"] @@ -443,22 +384,26 @@ uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" [[deps.Unicode]] uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" +[[deps.VectorizedKmers]] +git-tree-sha1 = "eaf497f1b720ae17dc52463dfe5fed5f74b0e4df" +uuid = "2ef45bd6-4c2b-43d8-88b3-40597287359a" +version = "0.8.1" +weakdeps = ["BioSequences"] + + [deps.VectorizedKmers.extensions] + BioSequencesExt = "BioSequences" + [[deps.Zlib_jll]] deps = ["Libdl"] uuid = "83775a58-1f1d-513f-b197-d71354ab007a" -version = "1.2.13+0" - -[[deps.libblastrampoline_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" -version = "5.8.0+0" +version = "1.2.13+1" [[deps.nghttp2_jll]] deps = ["Artifacts", "Libdl"] uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d" -version = "1.48.0+0" +version = "1.52.0+1" [[deps.p7zip_jll]] deps = ["Artifacts", "Libdl"] uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0" -version = "17.4.0+0" +version = "17.4.0+2" diff --git a/docs/Project.toml b/docs/Project.toml index 8697a28..a5e05f9 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,4 +1,5 @@ [deps] +BioMarkovChains = "f861b655-cb5f-42ce-b66a-341b542d4f2c" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" GeneFinder = "2bc6ee46-2158-4f5a-a720-91cb2d7cee64" LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589" diff --git a/docs/make.jl b/docs/make.jl index a2a0b9c..8ed9b1f 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -10,7 +10,7 @@ makedocs(; sitename = "GeneFinder.jl", format = Documenter.HTML(; mathengine=MathJax3(), - prettyurls = get(ENV, "CI", "false") == "true", + # prettyurls = get(ENV, "CI", "false") == "true", canonical = "https://camilogarciabotero.github.io/GeneFinder.jl", repolink = "https://github.com/camilogarciabotero/GeneFinder.jl", edit_link = "main", @@ -18,6 +18,8 @@ makedocs(; pages = [ "Home" => "index.md", "Finding ORFs" => "simplefinder.md", + "Wrtiting ORFs in files" => "iodocs.md", + # "A Simple coding rule" => "simplecodingrule.md", "Roadmap" => "roadmap.md", "API" => "api.md", ], diff --git a/docs/src/_index.md b/docs/src/_index.md deleted file mode 100644 index deab8be..0000000 --- a/docs/src/_index.md +++ /dev/null @@ -1,87 +0,0 @@ - -

-
A Gene Finder -framework for Julia.

- - - - - - - - - - -

- ------------------------------------------------------------------------- - -## Overview - -> This is a species-agnostic, algorithm extensible, sequence-anonymous -> (genome, metagenomes) *gene finder* library for the Julia Language. - -The main goal is to create a versatile module that enables apply -different implemented algorithm to DNA sequences. See, for instance, -BioAlignment implementations of different sequence alignment algorithms -(local, global, edit-distance). - -## Installation - -You can install GeneFinder from the julia REPL. Press `]` to enter pkg -mode, and enter the following: - - add GeneFinder - -If you are interested in the cutting edge of the development, please -check out the master branch to try new features before release. - -## Algorithms - -### Coding genes (CDS - ORFs) - -- ☒ [Simple - finder](https://camilogarciabotero.github.io/GeneFinder.jl/dev/simplefinder/) -- ☐ EasyGene -- ☐ GLIMMER -- ☐ Prodigal - Pyrodigal -- ☐ PHANOTATE -- ☐ k-mer based gene finders (?) -- ☐ Augustus (?) - -### Non-coding genes (RNA) - -- ☐ Infernal -- ☐ tRNAscan - -## Other features - -- ☐ parallelism SIMD ? -- ☐ memory management (?) -- ☐ specialized types - - ☒ Gene - - ☒ ORF - - ☒ CDS - - ☐ EukaryoticGene (?) - - ☐ ProkaryoticGene (?) - - ☐ Codon - - ☐ Intron - - ☐ Exon - - ☐ GFF –\> See other packages - - ☐ FASTX –\> See I/O in other packages - -## Compatibilities - -Must interact with or extend: - -- GenomicAnnotations.jl -- BioSequences.jl -- SequenceVariation.jl -- GenomicFeatures.jl -- FASTX.jl -- Kmers.jl - -## Contributing - -## Citing - diff --git a/docs/src/index.md b/docs/src/index.md index aa26b0e..8d1d39e 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -32,6 +32,9 @@ Downloads + + Aqua QA + @@ -63,7 +66,7 @@ add GeneFinder author = {Camilo García}, title = {GeneFinder.jl}, url = {https://github.com/camilogarciabotero/GeneFinder.jl}, - version = {v0.0.23}, + version = {v0.2.0}, year = {2022}, month = {11} } diff --git a/docs/src/iodocs.md b/docs/src/iodocs.md new file mode 100644 index 0000000..f698e35 --- /dev/null +++ b/docs/src/iodocs.md @@ -0,0 +1,96 @@ +## Writting ORF information into bioinformatic formats + +This package facilitates the creation of `FASTA`, `BED`, and `GFF` files, specifically extracting Open Reading Frame (ORF) information from `BioSequence` instances, particularly those of type `NucleicSeqOrView{A} where A`, and then writing the information into the desired format. + +Functionality: + +The package provides four distinct functions for writing files in different formats: + +| Function | Description | +|-------------------|--------------------------------------------------------| +| `write_orfs_fna` | Writes nucleotide sequences in FASTA format. | +| `write_orfs_faa` | Writes amino acid sequences in FASTA format. | +| `write_orfs_bed` | Outputs information in BED format. | +| `write_orfs_gff` | Generates files in GFF format. | + + +All these functions support processing both `BioSequence` instances and external `FASTA` files. In the case of a `BioSequence` instace into external files, simply provide the path to the `FASTA` file using a `String` to the path. To demonstrate the use of the `write_*` methods with a `BioSequence`, consider the following example: + +```julia +using BioSequences, GeneFinder + +# > 180195.SAMN03785337.LFLS01000089 -> finds only 1 gene in Prodigal (from Pyrodigal tests) +seq = dna"AACCAGGGCAATATCAGTACCGCGGGCAATGCAACCCTGACTGCCGGCGGTAACCTGAACAGCACTGGCAATCTGACTGTGGGCGGTGTTACCAACGGCACTGCTACTACTGGCAACATCGCACTGACCGGTAACAATGCGCTGAGCGGTCCGGTCAATCTGAATGCGTCGAATGGCACGGTGACCTTGAACACGACCGGCAATACCACGCTCGGTAACGTGACGGCACAAGGCAATGTGACGACCAATGTGTCCAACGGCAGTCTGACGGTTACCGGCAATACGACAGGTGCCAACACCAACCTCAGTGCCAGCGGCAACCTGACCGTGGGTAACCAGGGCAATATCAGTACCGCAGGCAATGCAACCCTGACGGCCGGCGACAACCTGACGAGCACTGGCAATCTGACTGTGGGCGGCGTCACCAACGGCACGGCCACCACCGGCAACATCGCGCTGACCGGTAACAATGCACTGGCTGGTCCTGTCAATCTGAACGCGCCGAACGGCACCGTGACCCTGAACACAACCGGCAATACCACGCTGGGTAATGTCACCGCACAAGGCAATGTGACGACTAATGTGTCCAACGGCAGCCTGACAGTCGCTGGCAATACCACAGGTGCCAACACCAACCTGAGTGCCAGCGGCAATCTGACCGTGGGCAACCAGGGCAATATCAGTACCGCGGGCAATGCAACCCTGACTGCCGGCGGTAACCTGAGC" +``` + +Once a `BioSequence` object has been instantiated, the `write_orfs_fna` function proves useful for generating a `FASTA` file containing the nucleotide sequences of the ORFs. Notably, the `write_orfs*` methods support either an `IOStream` or an `IOBuffer` as an output argument, allowing flexibility in directing the output either to a file or a buffer. In the following example, we demonstrate writing the output directly to a file. + +```julia +outfile = "LFLS01000089.fna" + +open(outfile, "w") do io + write_orfs_fna(seq, io) +end +``` + +```bash +cat LFLS01000089.fna + +>ORF01 id=01 start=29 stop=40 strand=+ frame=2 +ATGCAACCCTGA +>ORF02 id=02 start=137 stop=145 strand=+ frame=2 +ATGCGCTGA +>ORF03 id=03 start=164 stop=184 strand=+ frame=2 +ATGCGTCGAATGGCACGGTGA +>ORF04 id=04 start=173 stop=184 strand=+ frame=2 +ATGGCACGGTGA +>ORF05 id=05 start=236 stop=241 strand=+ frame=2 +ATGTGA +>ORF06 id=06 start=248 stop=268 strand=+ frame=2 +ATGTGTCCAACGGCAGTCTGA +>ORF07 id=07 start=362 stop=373 strand=+ frame=2 +ATGCAACCCTGA +>ORF08 id=08 start=470 stop=496 strand=+ frame=2 +ATGCACTGGCTGGTCCTGTCAATCTGA +>ORF09 id=09 start=551 stop=574 strand=+ frame=2 +ATGTCACCGCACAAGGCAATGTGA +>ORF10 id=10 start=569 stop=574 strand=+ frame=2 +ATGTGA +>ORF11 id=11 start=581 stop=601 strand=+ frame=2 +ATGTGTCCAACGGCAGCCTGA +>ORF12 id=12 start=695 stop=706 strand=+ frame=2 +ATGCAACCCTGA +``` + +### Combining `FASTX` for reading and writing fastas + +We can now combine the `FASTX` package with the function `write_orfs_faa` to write a `FASTA` file with the protein sequences of the translated ORFs obtained from an external `FASTA` file. + +```julia +infile = "test/data/NC_001884.fasta" +outfile = "test/data/NC_001884-orfs.faa" + +open(inputfile) do io + write_orfs_faa(infile, outfile) +end +``` + +```bash +head test/data/NC_001884-orfs.faa + +>ORF0001 id=0001 start=41 stop=145 strand=- frame=2 +MTQKRKGPIPAQFEITPILRFNFIFDLTATNSFH* +>ORF0002 id=0002 start=41 stop=172 strand=- frame=2 +MVLKDVIVNMTQKRKGPIPAQFEITPILRFNFIFDLTATNSFH* +>ORF0003 id=0003 start=41 stop=454 strand=- frame=2 +MSEHLSQKEKELKNKENFIFDKYESGIYSDELFLKRKAALDEEFKELQNAKNELNGLQDTQSEIDSNTVRNNINKIIDQYHIESSSEKKNELLRMVLKDVIVNMTQKRKGPIPAQFEITPILRFNFIFDL +TATNSFH* +>ORF0004 id=0004 start=41 stop=472 strand=- frame=2 +MKTKKQMSEHLSQKEKELKNKENFIFDKYESGIYSDELFLKRKAALDEEFKELQNAKNELNGLQDTQSEIDSNTVRNNINKIIDQYHIESSSEKKNELLRMVLKDVIVNMTQKRKGPIPAQFEITPILRF +NFIFDLTATNSFH* +>ORF0005 id=0005 start=41 stop=505 strand=- frame=2 +MLSKYEDDNSNMKTKKQMSEHLSQKEKELKNKENFIFDKYESGIYSDELFLKRKAALDEEFKELQNAKNELNGLQDTQSEIDSNTVRNNINKIIDQYHIESSSEKKNELLRMVLKDVIVNMTQKRKGPIP +AQFEITPILRFNFIFDLTATNSFH* +``` + +This could also be done to writting a `FASTA` file with the nucleotide sequences of the ORFs using the `write_orfs_fna` function. Similarly for the `BED` and `GFF` files using the `write_orfs_bed` and `write_orfs_gff` functions respectively. \ No newline at end of file diff --git a/docs/src/markovchains.html b/docs/src/markovchains.html deleted file mode 100644 index b451516..0000000 --- a/docs/src/markovchains.html +++ /dev/null @@ -1,445 +0,0 @@ - - - - - - - - - -markovchains - - - - - - - - - - - - - - - - - - - - - -
- -
- -
-
-

Representing DNA sequences as a Markov chain

-
- - - -
- - - - -
- - -
- -

One important step towards gene finding algorithms is to represent a DNA sequence as a Markov chain Figure 1. In this representation a DNA sequence of a reduced alphabet \(\scrA = {A, C, G, T}\) is draw as a four-vertex graph, where each letter of the \(\scrA\) is a state and the edges of the graph represent transitions from one nucleotide to another in a sequence. This is also consdired more specifically as a Discrete Markov chain @axelson-fisk_comparative_2015.

-
-
-

-
Figure 1: DNA sequence as a Markov chain
-
-
-

Now, several packages (e.g. MarkovChainsHammer, Markovians, etc.) in the Julia ecosystem have been

- -
- - -
- - - - - \ No newline at end of file diff --git a/docs/src/simplecodingrule.md b/docs/src/simplecodingrule.md new file mode 100644 index 0000000..6175547 --- /dev/null +++ b/docs/src/simplecodingrule.md @@ -0,0 +1,29 @@ +## The *log-odds ratio* decision rule + +The sequence probability given a transition probability model (eq. 2) could be used as the source of a sequence classification based on a decision rule to classify whether a sequence correspond to a model or another. Now, imagine we got two DNA sequence transition models, a CDS model and a No-CDS model. The *log-odds ratio* decision rule could be establish as: + +``` math +\begin{align} +S(X) = \log \frac{{P_C(X_1=i_1, \ldots, X_T=i_T)}}{{P_N(X_1=i_1, \ldots, X_T=i_T)}} \begin{cases} > \eta & \Rightarrow \text{coding} \\ < \eta & \Rightarrow \text{noncoding} \end{cases} +\end{align} +``` + +Where the ``P_{C}`` is the probability of the sequence given a +CDS model, ``P_{N}`` is the probability of the sequence given a +No-CDS model, the decision rule is finally based on whether the ratio is +greater or lesser than a given threshold *η* of significance level. + +In the GeneFinder we have implemented this rule and a couple of basic +transition probability models of CDS and No-CDS of *E. coli* from +Axelson-Fisk (2015) work. To check whether a random sequence could be +coding based on these decision we use the predicate `iscoding` with the +`ECOLICDS` and `ECOLINOCDS` models: + +``` julia +using GeneFinder, BioSequences +randseq = get_orfs_dna(randdnaseq(99))[1] # this will retrieved a random coding ORF + +iscoding(randseq, ECOLICDS, ECOLINOCDS) +``` + + true \ No newline at end of file diff --git a/docs/src/simplefinder.html b/docs/src/simplefinder.html deleted file mode 100644 index 9ff4243..0000000 --- a/docs/src/simplefinder.html +++ /dev/null @@ -1,555 +0,0 @@ - - - - - - - - - -simplefinder - - - - - - - - - - - - - - - - - - - -
- -
- - - -
-

A simple algorithm

-

The first implemented function is simplefinder a very non-restrictive ORF finder function that will catch all ORFs in a dedicated structure. Note that this will catch random ORFs not necesarily genes since it has no ORFs size or overlapping condition contraints. Thus it might consider aa"M*" a posible encoding protein from the resulting ORFs.

-
-
using BioSequences, GeneFinder
-
-# > 180195.SAMN03785337.LFLS01000089 -> finds only 1 gene in Prodigal (from Pyrodigal tests)
-seq = dna"AACCAGGGCAATATCAGTACCGCGGGCAATGCAACCCTGACTGCCGGCGGTAACCTGAACAGCACTGGCAATCTGACTGTGGGCGGTGTTACCAACGGCACTGCTACTACTGGCAACATCGCACTGACCGGTAACAATGCGCTGAGCGGTCCGGTCAATCTGAATGCGTCGAATGGCACGGTGACCTTGAACACGACCGGCAATACCACGCTCGGTAACGTGACGGCACAAGGCAATGTGACGACCAATGTGTCCAACGGCAGTCTGACGGTTACCGGCAATACGACAGGTGCCAACACCAACCTCAGTGCCAGCGGCAACCTGACCGTGGGTAACCAGGGCAATATCAGTACCGCAGGCAATGCAACCCTGACGGCCGGCGACAACCTGACGAGCACTGGCAATCTGACTGTGGGCGGCGTCACCAACGGCACGGCCACCACCGGCAACATCGCGCTGACCGGTAACAATGCACTGGCTGGTCCTGTCAATCTGAACGCGCCGAACGGCACCGTGACCCTGAACACAACCGGCAATACCACGCTGGGTAATGTCACCGCACAAGGCAATGTGACGACTAATGTGTCCAACGGCAGCCTGACAGTCGCTGGCAATACCACAGGTGCCAACACCAACCTGAGTGCCAGCGGCAATCTGACCGTGGGCAACCAGGGCAATATCAGTACCGCGGGCAATGCAACCCTGACTGCCGGCGGTAACCTGAGC"
-
-
726nt DNA Sequence:
-AACCAGGGCAATATCAGTACCGCGGGCAATGCAACCCTG…GCGGGCAATGCAACCCTGACTGCCGGCGGTAACCTGAGC
-
-
-
-

Finding all ORFs

-
-
simplefind(seq)
-
-
12-element Vector{ORF}:
- ORF(29:40, '+')
- ORF(137:145, '+')
- ORF(164:184, '+')
- ORF(173:184, '+')
- ORF(236:241, '+')
- ORF(248:268, '+')
- ORF(362:373, '+')
- ORF(470:496, '+')
- ORF(551:574, '+')
- ORF(569:574, '+')
- ORF(581:601, '+')
- ORF(695:706, '+')
-
-
-

Two other functions (cdsgenerator and proteingenerator) pass the sequence to simplefinder take the ORFs and act as generators of the sequence, so this way the can be collected in the REPL as an standard output or writeen into a file more conviniently using the FASTX IO system:

-
-
-

Generting cds and proteins with its ORF

-
-
cds = simplecdsgenerator(seq)
-
-
Base.Generator{Base.Iterators.Filter{GeneFinder.var"#4#6"{Int64}, Vector{ORF}}, GeneFinder.var"#3#5"{LongSequence{DNAAlphabet{4}}, LongSequence{DNAAlphabet{4}}}}(GeneFinder.var"#3#5"{LongSequence{DNAAlphabet{4}}, LongSequence{DNAAlphabet{4}}}(AACCAGGGCAATATCAGTACCGCGGGCAATGCAACCCTG…GCGGGCAATGCAACCCTGACTGCCGGCGGTAACCTGAGC, GCTCAGGTTACCGCCGGCAGTCAGGGTTGCATTGCCCGC…CAGGGTTGCATTGCCCGCGGTACTGATATTGCCCTGGTT), Base.Iterators.Filter{GeneFinder.var"#4#6"{Int64}, Vector{ORF}}(GeneFinder.var"#4#6"{Int64}(6), ORF[ORF(29:40, '+'), ORF(137:145, '+'), ORF(164:184, '+'), ORF(173:184, '+'), ORF(236:241, '+'), ORF(248:268, '+'), ORF(362:373, '+'), ORF(470:496, '+'), ORF(551:574, '+'), ORF(569:574, '+'), ORF(581:601, '+'), ORF(695:706, '+')]))
-
-
[i.sequence for i in cds]
-
-
12-element Vector{LongSequence{DNAAlphabet{4}}}:
- ATGCAACCCTGA
- ATGCGCTGA
- ATGCGTCGAATGGCACGGTGA
- ATGGCACGGTGA
- ATGTGA
- ATGTGTCCAACGGCAGTCTGA
- ATGCAACCCTGA
- ATGCACTGGCTGGTCCTGTCAATCTGA
- ATGTCACCGCACAAGGCAATGTGA
- ATGTGA
- ATGTGTCCAACGGCAGCCTGA
- ATGCAACCCTGA
-
-
-
-
protein = simpleprot_generator(seq)
-
-
Base.Generator{Base.Iterators.Filter{GeneFinder.var"#14#16"{Int64}, Vector{ORF}}, GeneFinder.var"#13#15"{Bool, BioSequences.GeneticCode, LongSequence{DNAAlphabet{4}}, LongSequence{DNAAlphabet{4}}}}(GeneFinder.var"#13#15"{Bool, BioSequences.GeneticCode, LongSequence{DNAAlphabet{4}}, LongSequence{DNAAlphabet{4}}}(false, The Standard Code, AACCAGGGCAATATCAGTACCGCGGGCAATGCAACCCTG…GCGGGCAATGCAACCCTGACTGCCGGCGGTAACCTGAGC, GCTCAGGTTACCGCCGGCAGTCAGGGTTGCATTGCCCGC…CAGGGTTGCATTGCCCGCGGTACTGATATTGCCCTGGTT), Base.Iterators.Filter{GeneFinder.var"#14#16"{Int64}, Vector{ORF}}(GeneFinder.var"#14#16"{Int64}(6), ORF[ORF(29:40, '+'), ORF(137:145, '+'), ORF(164:184, '+'), ORF(173:184, '+'), ORF(236:241, '+'), ORF(248:268, '+'), ORF(362:373, '+'), ORF(470:496, '+'), ORF(551:574, '+'), ORF(569:574, '+'), ORF(581:601, '+'), ORF(695:706, '+')]))
-
-
[i.sequence for i in protein]
-
-
12-element Vector{LongAA}:
- MQP*
- MR*
- MRRMAR*
- MAR*
- M*
- MCPTAV*
- MQP*
- MHWLVLSI*
- MSPHKAM*
- M*
- MCPTAA*
- MQP*
-
-
-
-
-

Combining FASTX to read a fasta record

-
-
using FASTX
-
-filename = "../../test/data/KK037166.fna"
-rdr = FASTA.Reader(open(filename))
-record = first(rdr)
-seq = sequence(record)
-dnaseq = LongDNA{4}(seq)
-[i.sequence for i in simpleprot_generator(dnaseq)]
-
-
-
-

Writting cds and proteins fastas

-
write_cds("cds.fasta", seq)
-
cat cds.fasta
->locus=29:40 strand=+
-ATGCAACCCTGA
->locus=137:145 strand=+
-ATGCGCTGA
->locus=164:184 strand=+
-ATGCGTCGAATGGCACGGTGA
->locus=173:184 strand=+
-ATGGCACGGTGA
->locus=236:241 strand=+
-ATGTGA
->locus=248:268 strand=+
-ATGTGTCCAACGGCAGTCTGA
->locus=362:373 strand=+
-ATGCAACCCTGA
->locus=470:496 strand=+
-ATGCACTGGCTGGTCCTGTCAATCTGA
->locus=551:574 strand=+
-ATGTCACCGCACAAGGCAATGTGA
->locus=569:574 strand=+
-ATGTGA
->locus=581:601 strand=+
-ATGTGTCCAACGGCAGCCTGA
->locus=695:706 strand=+
-ATGCAACCCTGA
-
write_proteins("proteins.fasta", seq)
-
cat proteins.fasta
->locus=29:40 strand=+
-MQP*
->locus=137:145 strand=+
-MR*
->locus=164:184 strand=+
-MRRMAR*
->locus=173:184 strand=+
-MAR*
->locus=236:241 strand=+
-M*
->locus=248:268 strand=+
-MCPTAV*
->locus=362:373 strand=+
-MQP*
->locus=470:496 strand=+
-MHWLVLSI*
->locus=551:574 strand=+
-MSPHKAM*
->locus=569:574 strand=+
-M*
->locus=581:601 strand=+
-MCPTAA*
->locus=695:706 strand=+
-MQP*
-
-
- -
- - -
- - - - \ No newline at end of file diff --git a/docs/src/simplefinder.md b/docs/src/simplefinder.md index 3a69ab3..faac22a 100644 --- a/docs/src/simplefinder.md +++ b/docs/src/simplefinder.md @@ -1,5 +1,5 @@ -## Finding complete and internal (overlapped) ORFs +## Finding complete and overlapped ORFs The first implemented function is `findorfs` a very non-restrictive ORF finder function that will catch all ORFs in a dedicated structure. Note that this will catch random ORFs not necesarily genes since it has no ORFs size or overlapping condition contraints. Thus it might consider `aa"M*"` a posible encoding protein from the resulting ORFs. @@ -29,10 +29,10 @@ findorfs(seq) ORF(695:706, '+', 2) ``` -Two other functions (`getorfdna` and `getorfaa`) pass the sequence to `findorfs` take the ORFs and act as generators of the sequence, so this way the can be `collect`ed in the REPL as an standard output or writteen into a file more conviniently using the `FASTX` IO system: +Two other functions (`get_orfs_dna` and `get_orfs_aa`) are implemented to get the ORFs in DNA and amino acid sequences, respectively. They use the `findorfs` function to first get the ORFs and then get the correspondance array of `BioSequence` objects. ```julia -getorfdna(seq) +get_orfs_dna(seq) 12-element Vector{LongSubSeq{DNAAlphabet{4}}}: ATGCAACCCTGA @@ -50,7 +50,7 @@ getorfdna(seq) ``` ```julia -getorfaa(seq) +get_orfs_aa(seq) 12-element Vector{LongSubSeq{AminoAcidAlphabet}}: MQP* @@ -67,94 +67,30 @@ getorfaa(seq) MQP* ``` -### Writting cds, proteins fastas, bed and gffs whether from a `LongSeq` or from a external fasta file. +## The ORF type -```julia -write_cds("cds.fasta", seq) -``` - -```bash -cat cds.fasta - ->location=29:40 strand=+ frame=2 -ATGCAACCCTGA ->location=137:145 strand=+ frame=2 -ATGCGCTGA ->location=164:184 strand=+ frame=2 -ATGCGTCGAATGGCACGGTGA ->location=173:184 strand=+ frame=2 -ATGGCACGGTGA ->location=236:241 strand=+ frame=2 -ATGTGA ->location=248:268 strand=+ frame=2 -ATGTGTCCAACGGCAGTCTGA ->location=362:373 strand=+ frame=2 -ATGCAACCCTGA ->location=470:496 strand=+ frame=2 -ATGCACTGGCTGGTCCTGTCAATCTGA ->location=551:574 strand=+ frame=2 -ATGTCACCGCACAAGGCAATGTGA ->location=569:574 strand=+ frame=2 -ATGTGA ->location=581:601 strand=+ frame=2 -ATGTGTCCAACGGCAGCCTGA ->location=695:706 strand=+ frame=2 -ATGCAACCCTGA -``` - -### Combining `FASTX` for reading and writing fastas +For convenience, the `ORF` type is more stringent in preventing the creation of incompatible instances. As a result, attempting to create an instance with incompatible parameters will result in an error. For instance, the following code snippet will trigger an error: ```julia -using FASTX - -write_proteins("test/data/NC_001884.fasta", "proteins.fasta") -``` - -```bash -head proteins.fasta - ->location=41:145 strand=- frame=2 -MTQKRKGPIPAQFEITPILRFNFIFDLTATNSFH* ->location=41:172 strand=- frame=2 -MVLKDVIVNMTQKRKGPIPAQFEITPILRFNFIFDLTATNSFH* ->location=41:454 strand=- frame=2 -MSEHLSQKEKELKNKENFIFDKYESGIYSDELFLKRKAALDEEFKELQNAKNELNGLQDTQSEIDSNTVRNNINKIIDQYHIESSSEKKNELLRMVLKDVIVNMTQKRKGPIPAQFEITPILRFNFIFDLTATNSFH* ->location=41:472 strand=- frame=2 -MKTKKQMSEHLSQKEKELKNKENFIFDKYESGIYSDELFLKRKAALDEEFKELQNAKNELNGLQDTQSEIDSNTVRNNINKIIDQYHIESSSEKKNELLRMVLKDVIVNMTQKRKGPIPAQFEITPILRFNFIFDLTATNSFH* ->location=41:505 strand=- frame=2 -MLSKYEDDNSNMKTKKQMSEHLSQKEKELKNKENFIFDKYESGIYSDELFLKRKAALDEEFKELQNAKNELNGLQDTQSEIDSNTVRNNINKIIDQYHIESSSEKKNELLRMVLKDVIVNMTQKRKGPIPAQFEITPILRFNFIFDLTATNSFH* -``` - -## The *log-odds ratio* decision rule - -The sequence probability given a transition probability model (eq. 2) -could be used as the source of a sequence classification based on a -decision rule to classify whether a sequence correspond to a model or -another. Now, imagine we got two DNA sequence transition models, a CDS -model and a No-CDS model. The *log-odds ratio* decision rule could be -establish as: - -``` math -\begin{align} -S(X) = \log \frac{{P_C(X_1=i_1, \ldots, X_T=i_T)}}{{P_N(X_1=i_1, \ldots, X_T=i_T)}} \begin{cases} > \eta & \Rightarrow \text{coding} \\ < \eta & \Rightarrow \text{noncoding} \end{cases} -\end{align} +ORF(1:10, '+', 4) + +ERROR: AssertionError: Invalid frame value. Frame must be 1, 2, or 3. +Stacktrace: + [1] ORF(location::UnitRange{Int64}, strand::Char, frame::Int64) + @ GeneFinder ~/.julia/dev/GeneFinder/src/types.jl:47 + [2] top-level scope + @ REPL[25]:1 ``` + +Similar behavior will be encountered when the strand is neither `+` nor `-`. This precautionary measure helps prevent the creation of invalid ORFs, ensuring greater stability and enabling the extension of its interface. For example, after creating a specific `ORF`, users can seamlessly iterate over a sequence of interest and verify whether the ORF is contained within the sequence. -Where the ``P_{C}`` is the probability of the sequence given a -CDS model, ``P_{N}`` is the probability of the sequence given a -No-CDS model, the decision rule is finally based on whether the ratio is -greater or lesser than a given threshold *η* of significance level. - -In the GeneFinder we have implemented this rule and a couple of basic -transition probability models of CDS and No-CDS of *E. coli* from -Axelson-Fisk (2015) work. To check whether a random sequence could be -coding based on these decision we use the predicate `iscoding` with the -`ECOLICDS` and `ECOLINOCDS` models: - -``` julia -randseq = getorfdna(randdnaseq(99))[1] # this will retrieved a random coding ORF +```julia +orf = ORF(137:145, '+', 2) +seq[orf] -iscoding(randseq, ECOLICDS, ECOLINOCDS) +9nt DNA Sequence: +ATGCGCTGA ``` - true \ No newline at end of file +!!! warning + It is still possible to create an `ORF` and pass it to a sequence that does not necessarily contain an actual open reading frame. This will be addressed in future versions of the package. But the benefit of having it is that it will retrieve the corresponding subsequence of the sequence in a convinient way (5' to 3') regardless of the strand. \ No newline at end of file diff --git a/docs/src/simplefinder.qmd b/docs/src/simplefinder.qmd index 32d9da7..44feb26 100644 --- a/docs/src/simplefinder.qmd +++ b/docs/src/simplefinder.qmd @@ -47,7 +47,7 @@ getproteins(seq) #| eval: false using FASTX -write_proteins("../../test/data/NC_001884.fasta", "proteins.fasta") +write_orfs_faa("../../test/data/NC_001884.fasta", "proteins.fasta") ``` ```bash @@ -101,7 +101,7 @@ ATGCAACCCTGA ``` ```julia -write_proteins("proteins.fasta", seq) +write_orfs_faa("proteins.fasta", seq) ``` ```bash diff --git a/references.bib b/references.bib deleted file mode 100644 index ab7c59b..0000000 --- a/references.bib +++ /dev/null @@ -1,15 +0,0 @@ -@book{axelson-fisk_comparative_2015, - address = {London}, - series = {Computational {Biology}}, - title = {Comparative {Gene} {Finding}}, - volume = {20}, - isbn = {978-1-4471-6692-4 978-1-4471-6693-1}, - url = {http://link.springer.com/10.1007/978-1-4471-6693-1}, - language = {en}, - urldate = {2022-12-18}, - publisher = {Springer London}, - author = {Axelson-Fisk, Marina}, - year = {2015}, - doi = {10.1007/978-1-4471-6693-1}, - keywords = {bioinformatics, programming, gene finding, software} -} \ No newline at end of file diff --git a/src/GeneFinder.jl b/src/GeneFinder.jl index 828e0b5..c9f0868 100644 --- a/src/GeneFinder.jl +++ b/src/GeneFinder.jl @@ -6,35 +6,36 @@ using BioSequences: NucleicAcidAlphabet, DNAAlphabet, + RNAAlphabet, AminoAcidAlphabet, LongDNA, LongAA, LongSequence, LongSubSeq, + NucleicSeqOrView, @biore_str, @dna_str, GeneticCode, reverse_complement, - randdnaseq + randdnaseq, + ncbi_trans_table, + translate - -using BioMarkovChains: BioMarkovChain -using FASTX: FASTA, sequence, FASTAReader +using FASTX: FASTAReader, FASTARecord, description, sequence using IterTools: takewhile, iterated using PrecompileTools: @setup_workload, @compile_workload -using StatsBase: countmap include("types.jl") export ORF include("algorithms/findorfs.jl") -export locationiterator, findorfs, getorfdna, getorfaa +export locationiterator, findorfs, get_orfs_dna, get_orfs_aa, record_orfs_fna, record_orfs_faa include("io.jl") -export write_cds, write_proteins, write_bed, write_gff +export write_orfs_fna, write_orfs_faa, write_orfs_bed, write_orfs_gff include("utils.jl") -export fasta_to_dna, nucleotidefreqs, hasprematurestop, dnaseqprobability, iscoding +export fasta_to_dna, hasprematurestop include("extended.jl") @@ -46,9 +47,9 @@ include("extended.jl") # all calls in this block will be precompiled, regardless of whether # they belong to your package or not (on Julia 1.8 and higher) findorfs(seq) - getorfdna(seq) - # getorfaa(seq) + get_orfs_dna(seq) + # get_orfs_aa(seq) end end -end +end \ No newline at end of file diff --git a/src/algorithms/findorfs.jl b/src/algorithms/findorfs.jl index c9e06f1..2ef8e8f 100644 --- a/src/algorithms/findorfs.jl +++ b/src/algorithms/findorfs.jl @@ -1,5 +1,6 @@ """ - locationiterator(sequence::LongSequence{DNAAlphabet{4}}; alternative_start::Bool=false) + locationiterator(sequence::NucleicSeqOrView{DNAAlphabet{N}}; alternative_start::Bool=false) where {N} + locationiterator(sequence::NucleicSeqOrView{RNAAlphabet{N}}; alternative_start::Bool=false) where {N} This is an iterator function that uses regular expressions to search the entire ORF (instead of start and stop codons) in a `LongSequence{DNAAlphabet{4}}` sequence. It uses an anonymous function that will find the first regularly expressed ORF. Then using this anonymous function it creates an iterator that will apply it until there is no other CDS. @@ -14,7 +15,10 @@ This is an iterator function that uses regular expressions to search the entire See more about the discussion [here](https://discourse.julialang.org/t/how-to-improve-a-generator-to-be-more-memory-efficient-when-it-is-collected/92932/8?u=camilogarciabotero) """ -function locationiterator(sequence::LongNucOrView{N}; alternative_start::Bool = false) where N +function locationiterator( + sequence::NucleicSeqOrView{DNAAlphabet{N}}; + alternative_start::Bool = false +) where {N} regorf = alternative_start ? biore"DTG(?:[N]{3})*?T(AG|AA|GA)"dna : biore"ATG(?:[N]{3})*?T(AG|AA|GA)"dna # regorf = alternative_start ? biore"DTG(?:[N]{3})*?T(AG|AA|GA)"dna : biore"ATG([N]{3})*T(AG|AA|GA)?"dna # an attempt to make it non PCRE non-determinsitic finder(x) = findfirst(regorf, sequence, first(x) + 1) # + 3 @@ -22,9 +26,10 @@ function locationiterator(sequence::LongNucOrView{N}; alternative_start::Bool = return itr end +const FRAMEDICT = Dict(0 => 3, 1 => 1, 2 => 2) # sets the int conversion for frame assignment later + """ - findorfs(sequence::LongSequence{DNAAlphabet{4}}; kwargs...)::Vector{ORF} - findorfs(sequence::String; kwargs...)::Vector{ORF} + findorfs(sequence::NucleicAlphabet{DNAAlphabet{N}}; alternative_start::Bool=false, min_len::Int64=6)::Vector{ORF} where {N} A simple implementation that finds ORFs in a DNA sequence. @@ -40,19 +45,21 @@ The `findorfs` function takes a LongSequence{DNAAlphabet{4}} sequence and return - `alternative_start::Bool=false`: If true will pass the extended start codons to search. This will increase 3x the exec. time. - `min_len::Int64=6`: Length of the allowed ORF. Default value allow `aa"M*"` a posible encoding protein from the resulting ORFs. """ -function findorfs(sequence::LongNucOrView{N}; alternative_start::Bool = false, min_len::Int64 = 6) where N - orfs = Vector{ORF}() +function findorfs( + sequence::NucleicSeqOrView{DNAAlphabet{N}}; + alternative_start::Bool = false, + min_len::Int64 = 6 +) where {N} + orfs = Vector{ORF}(undef, 0) reversedseq = reverse_complement(sequence) seqlen = length(sequence) - frames = Dict(0 => 3, 1 => 1, 2 => 2) - for strand in ('+', '-') seq = strand == '-' ? reversedseq : sequence - @inbounds for location in locationiterator(seq; alternative_start) + @inbounds for location in @views locationiterator(seq; alternative_start) if length(location) >= min_len - frame = strand == '+' ? frames[location.start % 3] : frames[(seqlen - location.stop + 1) % 3] + frame = strand == '+' ? FRAMEDICT[location.start % 3] : FRAMEDICT[(seqlen - location.stop + 1) % 3] push!(orfs, ORF(strand == '+' ? location : (seqlen - location.stop + 1):(seqlen - location.start + 1), strand, frame)) end end @@ -60,36 +67,16 @@ function findorfs(sequence::LongNucOrView{N}; alternative_start::Bool = false, m return sort(orfs) end -function findorfs(sequence::String; alternative_start::Bool = false, min_len::Int64 = 6) - sequence = LongSequence{DNAAlphabet{4}}(sequence) - orfs = Vector{ORF}() - reversedseq = reverse_complement(sequence) - seqlen = length(sequence) - - frames = Dict(0 => 3, 1 => 1, 2 => 2) - - for strand in ('+', '-') - seq = strand == '-' ? reversedseq : sequence - - @inbounds for location in locationiterator(seq; alternative_start) - if length(location) >= min_len - frame = strand == '+' ? frames[location.start % 3] : frames[(seqlen - location.stop + 1) % 3] - push!(orfs, ORF(strand == '+' ? location : (seqlen - location.stop + 1):(seqlen - location.start + 1), strand, frame)) - end - end - end - return sort(orfs) -end +#### get_orfs_* methods #### """ - getorfdna(input::LongSequence{DNAAlphabet{4}}, output::String; kwargs...) - getorfdna(input::String; kwargs...) ## for strings per se + get_orfs_dna(sequence::NucleicSeqOrView{DNAAlphabet{N}}, output::String; kwargs...) where {N} -This function takes a `LongSequence{DNAAlphabet{4}}` or `String` sequence and identifies the open reading frames (ORFs) using the `findorfs()` function. The function then extracts the DNA sequence of each ORF and stores it in a `Vector{LongSubSeq{DNAAlphabet{4}}}`. +This function takes a `NucleicSeqOrView{DNAAlphabet{N}}` sequence and identifies the open reading frames (ORFs) using the `findorfs()` function. The function then extracts the DNA sequence of each ORF and stores it in a `Vector{LongSubSeq{DNAAlphabet{4}}}`. # Arguments -- `input`: The input sequence as a `LongSequence{DNAAlphabet{4}}` or `String`. +- `sequence`: The input sequence as a `NucleicSeqOrView{DNAAlphabet{N}}` # Keyword Arguments @@ -97,35 +84,27 @@ This function takes a `LongSequence{DNAAlphabet{4}}` or `String` sequence and id - `min_len::Int64=6`: The minimum length of the allowed ORF. By default, it allows ORFs that can encode at least one amino acid (e.g., `aa"M*"`). """ -function getorfdna( - sequence::LongNucOrView{N}; +function get_orfs_dna( + sequence::NucleicSeqOrView{DNAAlphabet{N}}; alternative_start::Bool = false, min_len::Int64 = 6 -) where N +) where {N} orfs = findorfs(sequence; alternative_start, min_len) - seqs = Vector{LongSubSeq{DNAAlphabet{4}}}() - @inbounds for i in orfs - if i.strand == '+' - push!(seqs, @view sequence[i.location]) - else - newseq = reverse_complement(@view sequence[i.location]) - push!(seqs, newseq) - end + seqs = Vector{LongSubSeq{DNAAlphabet{N}}}(undef, length(orfs)) #Vector{}(undef, length(orfs)) # todo correct the output type + @inbounds for (i, orf) in enumerate(orfs) + seqs[i] = sequence[orf] end return seqs end -import BioSequences.standard_genetic_code - """ - getorfaa(input::LongSequence{DNAAlphabet{4}}; kwargs...) - getorfaa(input::String; kwargs...) ## for strings per se + get_orfs_aa(sequence::NucleicSeqOrView{DNAAlphabet{N}}; kwargs...) where {N} -This function takes a `LongSequence{DNAAlphabet{4}}` or `String` sequence and identifies the open reading frames (ORFs) using the `findorfs()` function. The function then translates each ORF into an amino acid sequence and stores it in a `Vector{LongSubSeq{AminoAcidAlphabet}}`. +This function takes a `NucleicSeqOrView{DNAAlphabet{N}}` sequence and identifies the open reading frames (ORFs) using the `findorfs()` function. The function then translates each ORF into an amino acid sequence and stores it in a `Vector{LongSubSeq{AminoAcidAlphabet}}`. # Arguments -- `input`: The input sequence as a `LongSequence{DNAAlphabet{4}}` or `String`. +- `sequence`: The input sequence as a `NucleicSeqOrView{DNAAlphabet{N}}` # Keyword Arguments @@ -133,21 +112,87 @@ This function takes a `LongSequence{DNAAlphabet{4}}` or `String` sequence and id - `min_len::Int64=6`: The minimum length of the allowed ORF. By default, it allows ORFs that can encode at least one amino acid (e.g., `aa"M*"`). """ -function getorfaa( - sequence::LongNucOrView{N}; +function get_orfs_aa( + sequence::NucleicSeqOrView{DNAAlphabet{N}}; alternative_start::Bool = false, - code::GeneticCode = standard_genetic_code, + code::GeneticCode = ncbi_trans_table[1], min_len::Int64 = 6 -) where N +) where {N} orfs = findorfs(sequence; alternative_start, min_len) - aas = Vector{LongSubSeq{AminoAcidAlphabet}}() - @inbounds for i in orfs - if i.strand == '+' - push!(aas, translate(@view sequence[i.location])) - else - newaa = translate(reverse_complement(@view sequence[i.location]); code) - push!(aas, newaa) - end + aas = Vector{LongSubSeq{AminoAcidAlphabet}}(undef, length(orfs)) + @inbounds for (i, orf) in enumerate(orfs) + aas[i] = translate(sequence[orf]; code) end return aas +end + +#### record_orfs_* methods #### + +""" + record_orfs_fna(sequence::NucleicSeqOrView{DNAAlphabet{N}}; kwargs...) where {N} + +Record Open Reading Frames (ORFs) in a nucleic acid sequence. + +# Arguments +- `sequence::NucleicSeqOrView{DNAAlphabet{N}}`: The nucleic acid sequence to search for ORFs. +- `alternative_start::Bool=false`: If `true`, consider alternative start codons for ORFs. +- `min_len::Int=6`: The minimum length of an ORF to be recorded. + +# Returns +An array of `FASTARecord` objects representing the identified ORFs. + +# Description +This function searches for Open Reading Frames (ORFs) in a given nucleic acid sequence. An ORF is a sequence of DNA that starts with a start codon and ends with a stop codon, without any other stop codons in between. By default, only the standard start codon (ATG) is considered, but if `alternative_start` is set to `true`, alternative start codons are also considered. The minimum length of an ORF to be recorded can be specified using the `min_len` argument. +""" +function record_orfs_fna( + sequence::NucleicSeqOrView{DNAAlphabet{N}}; + alternative_start::Bool = false, + min_len::Int64 = 6 +) where {N} + orfs = findorfs(sequence; alternative_start, min_len) + norfs = length(orfs) + padding = norfs < 10 ? length(string(norfs)) + 1 : length(string(norfs)) + records = FASTARecord[] + @inbounds for (index, orf) in enumerate(orfs) + id = string(lpad(string(index), padding, "0")) + header = "ORF$(id) id=$(id) start=$(orf.location.start) stop=$(orf.location.stop) strand=$(orf.strand) frame=$(orf.frame)" + record = FASTARecord(header, sequence[orf]) + push!(records, record) + end + return records +end + +""" + record_orfs_faa(sequence::NucleicSeqOrView{DNAAlphabet{N}}; kwargs...) where {N} + +This function takes a nucleic acid sequence and finds all open reading frames (ORFs) in the sequence. +An ORF is a sequence of DNA that starts with a start codon and ends with a stop codon. +The function returns a list of FASTA records, where each record represents an ORF in the sequence. + +# Arguments +- `sequence`: The nucleic acid sequence to search for ORFs. +- `alternative_start`: A boolean indicating whether to consider alternative start codons. Default is `false`. +- `code`: The genetic code to use for translation. Default is the NCBI translation table 1. +- `min_len`: The minimum length of an ORF. ORFs shorter than this length will be ignored. Default is 6. + +# Returns +- A list of FASTA records representing the ORFs found in the sequence. +""" +function record_orfs_faa( + sequence::NucleicSeqOrView{DNAAlphabet{N}}; + alternative_start::Bool = false, + code::GeneticCode = ncbi_trans_table[1], + min_len::Int64 = 6 +) where {N} + orfs = findorfs(sequence; alternative_start, min_len) + norfs = length(orfs) + padding = norfs < 10 ? length(string(norfs)) + 1 : length(string(norfs)) + records = FASTARecord[] + @inbounds for (index, orf) in enumerate(orfs) + id = string(lpad(string(index), padding, "0")) + header = "ORF$(id) id=$(id) start=$(orf.location.start) stop=$(orf.location.stop) strand=$(orf.strand) frame=$(orf.frame)" + record = FASTARecord(header, translate(sequence[orf]; code)) + push!(records, record) + end + return records end \ No newline at end of file diff --git a/src/extended.jl b/src/extended.jl index de38b2c..8fb8854 100644 --- a/src/extended.jl +++ b/src/extended.jl @@ -1,7 +1,8 @@ # Methods from main packages that expand their fuctions to this package structs -import Base: show, length, iterate, sort +import Base: length, iterate, sort, getindex Base.sort(v::Vector{<:ORF}) = sort(v, by = _orf_sort_key) +Base.getindex(sequence::NucleicSeqOrView{A}, orf::ORF) where {A} = orf.strand == '+' ? (@view sequence[orf.location]) : reverse_complement(@view sequence[orf.location]) function _orf_sort_key(orf::ORF) return (orf.location, orf.strand) diff --git a/src/io.jl b/src/io.jl index 57b9bdf..1b94522 100644 --- a/src/io.jl +++ b/src/io.jl @@ -1,65 +1,166 @@ """ - write_bed(input::LongSequence{DNAAlphabet{4}}, output::String; kwargs...) - write_bed(input::String, output::String; kwargs...) + write_orfs_bed(input::NucleicAcidAlphabet{DNAAlphabet{N}}, output::Union{IOStream, IOBuffer}; kwargs...) where {N} + write_orfs_bed(input::NucleicAcidAlphabet{DNAAlphabet{N}}, output::String; kwargs...) where {N} + write_orfs_bed(input::String, output::Union{IOStream, IOBuffer}; kwargs...) + write_orfs_bed(input::String, output::String; kwargs...) Write BED data to a file. -# Keywords - -- `alternative_start::Bool=false`: If true will pass the extended start codons to search. This will increase 3x the exec. time. -- `min_len::Int64=6`: Length of the allowed ORF. Default value allow `aa"M*"` a posible encoding protein from the resulting ORFs. +# Arguments +- `input::NucleicAcidAlphabet{DNAAlphabet{N}}`: The input DNA sequence. +- `output::String`: The output file path. +- `alternative_start::Bool=false`: If true, alternative start codons will be used when identifying CDSs. Default is `false`. +- `min_len::Int64=6`: The minimum length that a CDS must have in order to be included in the output file. Default is `6`. """ -function write_bed(input::LongSequence{DNAAlphabet{4}}, output::String; alternative_start = false, min_len = 6) +function write_orfs_bed( + input::NucleicSeqOrView{DNAAlphabet{N}}, + output::Union{IOStream, IOBuffer}; + alternative_start::Bool = false, + min_len::Int64 = 6 +) where {N} + orfs = findorfs(input; alternative_start, min_len) + @inbounds for orf in orfs + println(output, "$(orf.location.start)\t$(orf.location.stop)\t$(orf.strand)\t$(orf.frame)") + end +end + +function write_orfs_bed( + input::NucleicSeqOrView{DNAAlphabet{N}}, + output::String; + alternative_start::Bool = false, + min_len::Int64 = 6 +) where {N} + orfs = findorfs(input; alternative_start, min_len) open(output, "w") do f - @simd for i in findorfs(input; alternative_start, min_len) - write(f, "$(i.location.start)\t$(i.location.stop)\t$(i.strand)\t$(i.frame)\n") + @inbounds for orf in orfs + write(f, "$(orf.location.start)\t$(orf.location.stop)\t$(orf.strand)\t$(orf.frame)\n") end end end -function write_bed(input::String, output::String; alternative_start = false, min_len = 6) - dnaseq = fasta_to_dna(input)[1] +function write_orfs_bed( + input::String, + output::Union{IOStream, IOBuffer}; + alternative_start::Bool = false, + min_len::Int64 = 6 +) + input = fasta_to_dna(input)[1] # rewrite input to be a DNA sequence + orfs = findorfs(input; alternative_start, min_len) + @inbounds for orf in orfs + println(output, "$(orf.location.start)\t$(orf.location.stop)\t$(orf.strand)\t$(orf.frame)") + end +end +function write_orfs_bed( + input::String, + output::String; + alternative_start::Bool = false, + min_len::Int64 = 6 +) + input = fasta_to_dna(input)[1] + orfs = findorfs(input; alternative_start, min_len) open(output, "w") do f - @simd for i in findorfs(dnaseq; alternative_start, min_len) - write(f, "$(i.location.start)\t$(i.location.stop)\t$(i.strand)\t$(i.frame)\n") + @inbounds for orf in orfs + write(f, "$(orf.location.start)\t$(orf.location.stop)\t$(orf.strand)\t$(orf.frame)\n") end end end """ - write_cds(input::LongSequence{DNAAlphabet{4}}, output::String; kwargs...) - write_cds(input::String, output::String; kwargs...) + write_orfs_fna(input::NucleicSeqOrView{DNAAlphabet{N}}, output::Union{IOStream, IOBuffer}; kwargs...) where {N} + write_orfs_fna(input::NucleicSeqOrView{DNAAlphabet{N}}, output::String; kwargs...) where {N} + write_orfs_fna(input::String, output::Union{IOStream, IOBuffer}; kwargs...) + write_orfs_fna(input::String, output::String; kwargs...) Write a file containing the coding sequences (CDSs) of a given DNA sequence to the specified file. # Keywords -- `alternative_start`: A boolean value indicating whether alternative start codons should be used when identifying CDSs. Default is `false`. -- `min_len`: An integer representing the minimum length that a CDS must have in order to be included in the output file. Default is `6`. +- `alternative_start::Bool=false`: If true, alternative start codons will be used when identifying CDSs. Default is `false`. +- `min_len::Int64=6`: The minimum length that a CDS must have in order to be included in the output file. Default is `6`. +# Examples + +```julia +filename = "output.fna" + +seq = dna"ATGATGCATGCATGCATGCTAGTAACTAGCTAGCTAGCTAGTAA" + +open(filename, "w") do file + write_orfs_fna(seq, file) +end +``` """ -function write_cds(input::LongSequence{DNAAlphabet{4}}, output::String; alternative_start = false, min_len = 6) +function write_orfs_fna( + input::NucleicSeqOrView{DNAAlphabet{N}}, + output::Union{IOStream, IOBuffer}; + alternative_start::Bool = false, + min_len::Int64 = 6 +) where {N} + orfs = findorfs(input; alternative_start, min_len) + norfs = length(orfs) + padding = norfs < 10 ? length(string(norfs)) + 1 : length(string(norfs)) + @inbounds for (i, orf) in enumerate(orfs) + id = string(lpad(string(i), padding, "0")) + println(output, ">ORF$(id) id=$(id) start=$(orf.location.start) stop=$(orf.location.stop) strand=$(orf.strand) frame=$(orf.frame)\n$(input[orf])") # i.strand == '+' ? input[i.location] : reverse_complement(@view input[i.location]) + end +end + +function write_orfs_fna( + input::NucleicSeqOrView{DNAAlphabet{N}}, + output::String; + alternative_start::Bool = false, + min_len::Int64 = 6 +) where {N} + orfs = findorfs(input; alternative_start, min_len) + norfs = length(orfs) + padding = norfs < 10 ? length(string(norfs)) + 1 : length(string(norfs)) open(output, "w") do f - for i in findorfs(input; alternative_start, min_len) - # sequence = i.strand == '+' ? input[i.location] : reverse_complement(@view input[i.location]) - write(f, ">location=$(i.location) strand=$(i.strand) frame=$(i.frame)\n$(i.strand == '+' ? dnaseq[i.location] : reverse_complement(@view dnaseq[i.location]))\n") + for (i, orf) in enumerate(orfs) + id = string(lpad(string(i), padding, "0")) + write(f, ">ORF$(id) id=$(id) start=$(orf.location.start) stop=$(orf.location.stop) strand=$(orf.strand) frame=$(orf.frame)\n$(input[orf])\n") # i.strand == '+' ? input[i.location] : reverse_complement(@view input[i.location]) end end end -function write_cds(input::String, output::String; alternative_start = false, min_len = 6) - dnaseq = fasta_to_dna(input)[1] +function write_orfs_fna( + input::String, + output::Union{IOStream, IOBuffer}; + alternative_start::Bool = false, + min_len::Int64 = 6 +) + input = fasta_to_dna(input)[1] + orfs = findorfs(input; alternative_start, min_len) + norfs = length(orfs) + padding = norfs < 10 ? length(string(norfs)) + 1 : length(string(norfs)) + @inbounds for (i, orf) in enumerate(orfs) + id = string(lpad(string(i), padding, "0")) + println(output, ">ORF$(id) id=$(id) start=$(orf.location.start) stop=$(orf.location.stop) strand=$(orf.strand) frame=$(orf.frame)\n$(input[orf])") # i.strand == '+' ? input[i.location] : reverse_complement(@view input[i.location]) + end +end +function write_orfs_fna( + input::String, + output::String; + alternative_start::Bool = false, + min_len::Int64 = 6 +) + input = fasta_to_dna(input)[1] + orfs = findorfs(input; alternative_start, min_len) + norfs = length(orfs) + padding = norfs < 10 ? length(string(norfs)) + 1 : length(string(norfs)) open(output, "w") do f - for i in findorfs(dnaseq; alternative_start, min_len) - # sequence = i.strand == '+' ? dnaseq[i.location] : reverse_complement(@view dnaseq[i.location]) - write(f, ">location=$(i.location) strand=$(i.strand) frame=$(i.frame)\n$(i.strand == '+' ? dnaseq[i.location] : reverse_complement(@view dnaseq[i.location]))\n") + for (i, orf) in enumerate(orfs) + id = string(lpad(string(i), padding, "0")) + write(f, ">ORF$(id) id=$(id) start=$(orf.location.start) stop=$(orf.location.stop) strand=$(orf.strand) frame=$(orf.frame)\n$(input[orf])\n") # i.strand == '+' ? input[i.location] : reverse_complement(@view input[i.location]) end end end """ - write_proteins(input::LongSequence{DNAAlphabet{4}}, output::String; kwargs...) + write_orfs_faa(input::LongSequence{DNAAlphabet{4}}, output::Union{IOStream, IOBuffer}; kwargs...) where {N} + write_orfs_faa(input::LongSequence{DNAAlphabet{4}}, output::String; kwargs...) where {N} + write_orfs_faa(input::String, output::Union{IOStream, IOBuffer}; kwargs...) + write_orfs_faa(input::String, output::String; kwargs...) Write the protein sequences encoded by the coding sequences (CDSs) of a given DNA sequence to the specified file. @@ -68,75 +169,183 @@ Write the protein sequences encoded by the coding sequences (CDSs) of a given DN - `code::GeneticCode=BioSequences.standard_genetic_code`: The genetic code by which codons will be translated. See `BioSequences.ncbi_trans_table` for more info. - `alternative_start::Bool=false`: If true will pass the extended start codons to search. This will increase 3x the exec. time. - `min_len::Int64=6`: Length of the allowed ORF. Default value allow `aa"M*"` a posible encoding protein from the resulting ORFs. +# Examples + +```julia +filename = "output.faa" + +seq = dna"ATGATGCATGCATGCATGCTAGTAACTAGCTAGCTAGCTAGTAA" + +open(filename, "w") do file + write_orfs_faa(seq, file) +end +``` """ -function write_proteins( - input::LongSequence{DNAAlphabet{4}}, +function write_orfs_faa( + input::NucleicSeqOrView{DNAAlphabet{N}}, + output::Union{IOStream, IOBuffer}; + alternative_start::Bool = false, + code::GeneticCode = ncbi_trans_table[1], + min_len::Int64 = 6 +) where {N} + orfs = findorfs(input; alternative_start, min_len) + norfs = length(orfs) + padding = norfs < 10 ? length(string(norfs)) + 1 : length(string(norfs)) + @inbounds for (i, orf) in enumerate(orfs) + id = string(lpad(string(i), padding, "0")) + println(output, ">ORF$(id) id=$(id) start=$(orf.location.start) stop=$(orf.location.stop) strand=$(orf.strand) frame=$(orf.frame)\n$(translate(input[orf]; code))") + end +end + +function write_orfs_faa( + input::NucleicSeqOrView{DNAAlphabet{N}}, output::String; - alternative_start = false, - code::GeneticCode = BioSequences.standard_genetic_code, - min_len = 6, -) + alternative_start::Bool = false, + code::GeneticCode = ncbi_trans_table[1], + min_len::Int64 = 6, +) where {N} + orfs = findorfs(input; alternative_start, min_len) + norfs = length(orfs) + padding = norfs < 10 ? length(string(norfs)) + 1 : length(string(norfs)) open(output, "w") do f - for i in findorfs(input; alternative_start, code, min_len) - translation = i.strand == '+' ? translate(dnaseq[i.location]) : translate(reverse_complement(@view dnaseq[i.location])) - write(f, ">location=$(i.location) strand=$(i.strand) frame=$(i.frame)\n$(translation)\n") + for (i, orf) in enumerate(orfs) + id = string(lpad(string(i), padding, "0")) + write(f, ">ORF$(id) id=$(id) start=$(orf.location.start) stop=$(orf.location.stop) strand=$(orf.strand) frame=$(orf.frame)\n$(translate(input[orf]; code))\n") end end end -function write_proteins( +function write_orfs_faa( + input::String, + output::Union{IOStream, IOBuffer}; + alternative_start::Bool = false, + code::GeneticCode = ncbi_trans_table[1], + min_len::Int64 = 6 +) + input = fasta_to_dna(input)[1] + orfs = findorfs(input; alternative_start, min_len) + norfs = length(orfs) + padding = norfs < 10 ? length(string(norfs)) + 1 : length(string(norfs)) + @inbounds for (i, orf) in enumerate(orfs) + id = string(lpad(string(i), padding, "0")) + println(output, ">ORF$(id) id=$(id) start=$(orf.location.start) stop=$(orf.location.stop) strand=$(orf.strand) frame=$(orf.frame)\n$(translate(input[orf]; code))") + end +end + +function write_orfs_faa( input::String, output::String; - alternative_start = false, - min_len = 6, + alternative_start::Bool = false, + code::GeneticCode = ncbi_trans_table[1], + min_len::Int64 = 6, ) - dnaseq = fasta_to_dna(input)[1] - + input = fasta_to_dna(input)[1] + orfs = findorfs(input; alternative_start, min_len) + norfs = length(orfs) + padding = norfs < 10 ? length(string(norfs)) + 1 : length(string(norfs)) open(output, "w") do f - for i in findorfs(dnaseq; alternative_start, min_len) - translation = i.strand == '+' ? translate(dnaseq[i.location]) : translate(reverse_complement(@view dnaseq[i.location])) - write(f, ">location=$(i.location) strand=$(i.strand) frame=$(i.frame)\n$(translation)\n") + for (i, orf) in enumerate(orfs) + id = string(lpad(string(i), padding, "0")) + write(f, ">ORF$(id) id=$(id) start=$(orf.location.start) stop=$(orf.location.stop) strand=$(orf.strand) frame=$(orf.frame)\n$(translate(input[orf]; code))\n") end end end """ - write_gff(input::LongSequence{DNAAlphabet{4}}, output::String; kwargs...) - write_gff(input::String, output::String; kwargs...) + write_orfs_gff(input::NucleicSeqOrView{DNAAlphabet{N}}, output::Union{IOStream, IOBuffer}; kwargs...) where {N} + write_orfs_gff(input::NucleicSeqOrView{DNAAlphabet{N}}, output::String; kwargs...) where {N} + write_orfs_gff(input::String, output::Union{IOStream, IOBuffer}; kwargs...) + write_orfs_gff(input::String, output::String; kwargs...) Write GFF data to a file. -# Keywords +# Arguments +- `input`: The input DNA sequence. +- `output`: The output file to write the GFF data to. +- `alternative_start::Bool=false`: If true, extended start codons will be considered during the search, increasing the execution time. +- `min_len::Int64=6`: The minimum length of the allowed ORF. The default value allows for possible encoding proteins with the `aa"M*"` sequence. -- `alternative_start::Bool=false`: If true will pass the extended start codons to search. This will increase 3x the exec. time. -- `min_len::Int64=6`: Length of the allowed ORF. Default value allow `aa"M*"` a posible encoding protein from the resulting ORFs. """ -function write_gff(input::String, output::String; alternative_start = false, min_len = 6) - dnaseq = fasta_to_dna(input)[1] +function write_orfs_gff( + input::NucleicSeqOrView{A}, + output::Union{IOStream, IOBuffer}; + alternative_start::Bool = false, + min_len::Int64 = 6 +) where {A} + orfs = findorfs(input; alternative_start, min_len) + norfs = length(orfs) + padding = norfs < 10 ? length(string(norfs)) + 1 : length(string(norfs)) + println(output, "##gff-version 3\n##sequence-region Chr 1 $(length(input))") + for (i, orf) in enumerate(orfs) + id = string("ORF", lpad(string(i), padding, "0")) + println( + output, + "Chr\t.\tORF\t$(orf.location.start)\t$(orf.location.stop)\t.\t$(orf.strand)\t.\tID=$(id);Name=$(id);Frame=$(orf.frame)", + ) + end +end +function write_orfs_gff( + input::NucleicSeqOrView{A}, + output::String; + alternative_start::Bool = false, + min_len::Int64 = 6 +) where {A} + orfs = findorfs(input; alternative_start, min_len) + norfs = length(orfs) + padding = norfs < 10 ? length(string(norfs)) + 1 : length(string(norfs)) open(output, "w") do f - write(f, "##gff-version 3\n##sequence-region Chr 1 $(length(dnaseq))\n") - for (index, i) in enumerate(findorfs(dnaseq; alternative_start, min_len)) - id = string("ORF", lpad(string(index), 5, "0")) + write(f, "##gff-version 3\n##sequence-region Chr 1 $(length(input))\n") + for (i, orf) in enumerate(orfs) + id = string("ORF", lpad(string(i), padding, "0")) write( f, - "Chr\t.\tORF\t$(i.location.start)\t$(i.location.stop)\t.\t$(i.strand)\t.\tID=$(id);Name=$(id);Frame=$(i.frame)\n", + "Chr\t.\tORF\t$(orf.location.start)\t$(orf.location.stop)\t.\t$(orf.strand)\t.\tID=$(id);Name=$(id);Frame=$(orf.frame)\n", ) end end end -function write_gff(input::LongSequence{DNAAlphabet{4}}, output::String; alternative_start = false, min_len = 6) +function write_orfs_gff( + input::String, + output::Union{IOStream, IOBuffer}; + alternative_start::Bool = false, + min_len::Int64 = 6 +) + input = fasta_to_dna(input)[1] + orfs = findorfs(input; alternative_start, min_len) + norfs = length(orfs) + padding = norfs < 10 ? length(string(norfs)) + 1 : length(string(norfs)) + println(output, "##gff-version 3\n##sequence-region Chr 1 $(length(input))") + for (i, orf) in enumerate(orfs) + id = string("ORF", lpad(string(i), padding, "0")) + println( + output, + "Chr\t.\tORF\t$(orf.location.start)\t$(orf.location.stop)\t.\t$(orf.strand)\t.\tID=$(id);Name=$(id);Frame=$(orf.frame)", + ) + end +end + +function write_orfs_gff( + input::String, + output::String; + alternative_start::Bool = false, + min_len::Int64 = 6 +) + input = fasta_to_dna(input)[1] + orfs = findorfs(input; alternative_start, min_len) + norfs = length(orfs) + padding = norfs < 10 ? length(string(norfs)) + 1 : length(string(norfs)) open(output, "w") do f write(f, "##gff-version 3\n##sequence-region Chr 1 $(length(input))\n") - for (index, i) in enumerate(findorfs(input; alternative_start, min_len)) - id = string("ORF", lpad(string(index), 5, "0")) + for (i, orf) in enumerate(orfs) + id = string("ORF", lpad(string(i), padding, "0")) write( f, - "Chr\t.\tORF\t$(i.location.start)\t$(i.location.stop)\t.\t$(i.strand)\t.\tID=$(id);Name=$(id);Frame=$(i.frame)\n", + "Chr\t.\tORF\t$(orf.location.start)\t$(orf.location.stop)\t.\t$(orf.strand)\t.\tID=$(id);Name=$(id);Frame=$(orf.frame)\n", ) end end end -# Chr needs to be changed to fit well on IGV \ No newline at end of file +# Chr needs to be changed for fasta (?) or the name of seq (?) to fit well on IGV \ No newline at end of file diff --git a/src/types.jl b/src/types.jl index 7e00ac9..503af3a 100644 --- a/src/types.jl +++ b/src/types.jl @@ -1,8 +1,5 @@ # Structs associated with gene models - -abstract type Gene end - -const LongNucOrView{N} = Union{LongSequence{<:NucleicAcidAlphabet{N}},LongSubSeq{<:NucleicAcidAlphabet{N}}} +abstract type AbstractGene end """ struct GeneFeatures @@ -19,21 +16,21 @@ This is the main Gene struct, based on the fields that could be found in a GFF3, The idea is correct the frame and attributes that will have something like a possible list (id=Char;name=;locus_tag). The `write` and `get` functions should have a dedicated method for this struct. """ -struct GeneFeatures <: Gene - seqname::String - start::Int64 - stop::Int64 - score::Float64 - strand::Char - frame::Int8 # But maybe a Union to allow empty when reading a GFF? - attribute::String # Should be a Dict perhaps -end +# struct GeneFeatures <: AbstractGene +# seqname::String +# start::Int64 +# stop::Int64 +# score::Float64 +# strand::Char +# frame::Int8 # But maybe a Union to allow empty when reading a GFF? +# attribute::Dict # Should be a Dict perhaps +# end """ struct ORF location::UnitRange{Int64} strand::Char - frame::Integer + frame::Int end The ORF struct represents an open reading frame in a DNA sequence. It has two fields: @@ -41,58 +38,14 @@ The ORF struct represents an open reading frame in a DNA sequence. It has two fi - `location`: which is a UnitRange{Int64} indicating the start and end locations of the ORF in the sequence - `strand`: is a Char type indicating whether the ORF is on the forward ('+') or reverse ('-') strand of the sequence. """ -struct ORF <: Gene - location::UnitRange{Int64} # Note that it is also called position for gene struct in GenomicAnotations +struct ORF <: AbstractGene + location::UnitRange{Int64} # Note that it is also called position for gene struct in GenomicAnotations strand::Char - frame::Integer # 1, 2, or 3 -end - -##### The following implementation is from https://biojulia.dev/BioSequences.jl/stable/interfaces/ ##### -# """ -# Codon <: BioSequence{DNAAlphabet{2}} - -# A `Struct` representing a codon, which is a subtype of `BioSequence` with -# an `Alphabet` of type `DNAAlphabet{2}`. It has a single field `x` of type -# `UInt8`. This was implemente in The following implementation is from https://biojulia.dev/BioSequences.jl/stable/interfaces/ -# """ -# struct Codon <: BioSequence{DNAAlphabet{2}} -# x::UInt8 - -# function Codon(iterable) -# length(iterable) == 3 || error("Must have length 3") -# x = zero(UInt) -# for (i, nt) in enumerate(iterable) -# x |= BioSequences.encode(Alphabet(Codon), convert(DNA, nt)) << (6 - 2i) -# end -# new(x % UInt8) -# end -# end - -# Base.copy(seq::Codon) = Codon(seq.x) -# Base.count(codon::Codon, sequence::LongSequence{DNAAlphabet{4}}) = count(codon, sequence) -# Base.length(::Codon) = 3 + frame::Int # Use Int64 instead of Int -# function Base.count(codons::Vector{Codon}, sequence::LongSequence{DNAAlphabet{4}}) -# a = 0 -# @inbounds for i in codons -# a += count(i, sequence) -# end -# return a -# end - -# # has_interface(BioSequence, Codon, [DNA_C, DNA_T, DNA_G], false) - -# function extract_encoded_element(x::Codon, i::Int) -# ((x.x >>> (6 - 2i)) & 3) % UInt -# end - -# function translate( -# ntseq::Codon; -# code::GeneticCode = standard_genetic_code, -# allow_ambiguous_codons = true, -# alternative_start = false, -# ) -# translate(ntseq; code, allow_ambiguous_codons, alternative_start) -# end - -##### ---------------------------------------- ##### + function ORF(location::UnitRange{Int64}, strand::Char, frame::Int) + @assert frame in (1, 2, 3) "Invalid frame value. Frame must be 1, 2, or 3." + @assert strand in ('+', '-') "Invalid strand value. Strand must be '+' or '-'." + new(location, strand, frame) + end +end \ No newline at end of file diff --git a/src/utils.jl b/src/utils.jl index b2caa66..71ba839 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -16,37 +16,6 @@ end # end # end -""" - nucleotidefreqs(sequence::LongSequence{DNAAlphabet{4}}) -> Dict{DNA, Float64} - -Calculate the frequency of each nucleotide in a DNA sequence. - -# Arguments -- `sequence::LongSequence{DNAAlphabet{4}}`: A `LongSequence{DNAAlphabet{4}}` sequence. - -# Returns -A dictionary with each nucleotide in the sequence as a key, and its frequency as a value. - -# Example -``` -seq = dna"CCTCCCGGACCCTGGGCTCGGGAC" - -nucleotidefreqs(seq) - -Dict{DNA, Float64} with 4 entries: -DNA_T => 0.125 -DNA_A => 0.0833333 -DNA_G => 0.333333 -DNA_C => 0.458333 -``` -""" -function nucleotidefreqs(sequence::LongNucOrView{4})::Dict{DNA,Float64} - counts = countmap(sequence) - T = length(sequence) - F = Dict(i => counts[i] / T for i in keys(counts)) - return F -end - """ hasprematurestop(sequence::LongNucOrView{4})::Bool @@ -54,7 +23,7 @@ Determine whether the `sequence` of type `LongSequence{DNAAlphabet{4}}` contains Returns a boolean indicating whether the `sequence` has more than one stop codon. """ -function hasprematurestop(sequence::LongNucOrView{N})::Bool where N +function hasprematurestop(sequence::NucleicSeqOrView{A})::Bool where {A} stopcodons = [LongDNA{4}("TAA"), LongDNA{4}("TAG"), LongDNA{4}("TGA")] # Create a set of stop codons diff --git a/test/findorfstest.jl b/test/findorfstest.jl index 57c7607..6ff0bcc 100644 --- a/test/findorfstest.jl +++ b/test/findorfstest.jl @@ -39,20 +39,21 @@ @test length(NC_001416_orfs) == 885 end -@testset "getorfdna" begin +@testset "get_orfs_dna" begin seq01 = dna"ATGATGCATGCATGCATGCTAGTAACTAGCTAGCTAGCTAGTAA" - orfseqs = getorfdna(seq01) + sseq01 = "ATGATGCATGCATGCATGCTAGTAACTAGCTAGCTAGCTAGTAA" + orfseqs = get_orfs_dna(seq01) @test length(orfseqs) == 5 @test orfseqs[1] == dna"ATGATGCATGCATGCATGCTAGTAACTAGCTAG" end -# @testset "getorfaa" begin +@testset "get_orfs_aa" begin -# seq01 = dna"ATGATGCATGCATGCATGCTAGTAACTAGCTAGCTAGCTAGTAA" -# aas = getorfaa(seq01) + seq01 = dna"ATGATGCATGCATGCATGCTAGTAACTAGCTAGCTAGCTAGTAA" + aas = get_orfs_aa(seq01) -# @test length(aas) == 5 -# @test aas[1] == aa"MMHACMLVTS*" -# end \ No newline at end of file + @test length(aas) == 5 + @test aas[1] == aa"MMHACMLVTS*" +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index ffb199c..e9211a3 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,6 +7,21 @@ using GeneFinder include("findorfstest.jl") +using Aqua: + test_ambiguities, + test_persistent_tasks, + test_piracies, + test_stale_deps, + test_unbound_args, + test_undefined_exports + +test_ambiguities(GeneFinder) +test_persistent_tasks(GeneFinder) +test_piracies(GeneFinder) +test_stale_deps(GeneFinder) +test_unbound_args(GeneFinder) +test_undefined_exports(GeneFinder) + end # @testset "GeneFinder.jl" begin # # Write your tests here.