Moving ORF as an AbstractGenomicInterval{T}

Manifest.toml

@@ -2,14 +2,17 @@
 
 julia_version = "1.10.3"
 manifest_format = "2.0"
-project_hash = "b881dd48d80d6b19dfa1756a0b2d6e2e42d93e33"
+project_hash = "514456eec9673906fe896a3ae36685951a71c19f"
 
 [[deps.Automa]]
 deps = ["PrecompileTools", "TranscodingStreams"]
 git-tree-sha1 = "588e0d680ad1d7201d4c6a804dcb1cd9cba79fbb"
 uuid = "67c07d97-cdcb-5c2c-af73-a7f9c32a568b"
 version = "1.0.3"
 
+[[deps.Base64]]
+uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+
 [[deps.BioGenerics]]
 deps = ["TranscodingStreams"]
 git-tree-sha1 = "7bbc085aebc6faa615740b63756e4986c9e85a70"
@@ -42,6 +45,25 @@ git-tree-sha1 = "e32a61f028b823a172c75e26865637249bb30dff"
 uuid = "3c28c6f8-a34d-59c4-9654-267d177fcfa9"
 version = "5.1.3"
 
+[[deps.Compat]]
+deps = ["TOML", "UUIDs"]
+git-tree-sha1 = "b1c55339b7c6c350ee89f2c1604299660525b248"
+uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
+version = "4.15.0"
+
+    [deps.Compat.extensions]
+    CompatLinearAlgebraExt = "LinearAlgebra"
+
+    [deps.Compat.weakdeps]
+    Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
+    LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+
+[[deps.DataStructures]]
+deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
+git-tree-sha1 = "1d0a14036acb104d9e89698bd408f63ab58cdc82"
+uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+version = "0.18.20"
+
 [[deps.Dates]]
 deps = ["Printf"]
 uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
@@ -56,11 +78,35 @@ weakdeps = ["BioSequences"]
     [deps.FASTX.extensions]
     BioSequencesExt = "BioSequences"
 
+[[deps.GenomicFeatures]]
+deps = ["BioGenerics", "DataStructures", "IntervalTrees"]
+git-tree-sha1 = "720844f71f118dd2ab4898a89b2b4feca7730d81"
+uuid = "899a7d2d-5c61-547b-bef9-6698a8d05446"
+version = "3.0.0"
+
+[[deps.InteractiveUtils]]
+deps = ["Markdown"]
+uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+
+[[deps.IntervalTrees]]
+git-tree-sha1 = "dc3b97bb5c9cb7c437f74027309f2c2f09a82aaf"
+uuid = "524e6230-43b7-53ae-be76-1e9e4d08d11b"
+version = "1.1.0"
+
 [[deps.IterTools]]
 git-tree-sha1 = "42d5f897009e7ff2cf88db414a389e5ed1bdd023"
 uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
 version = "1.10.0"
 
+[[deps.Markdown]]
+deps = ["Base64"]
+uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
+
+[[deps.OrderedCollections]]
+git-tree-sha1 = "dfdf5519f235516220579f949664f1bf44e741c5"
+uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
+version = "1.6.3"
+
 [[deps.PrecompileTools]]
 deps = ["Preferences"]
 git-tree-sha1 = "5aa36f7049a63a1528fe8f7c3f2113413ffd4e1f"
@@ -117,6 +163,10 @@ 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"
 

Project.toml

@@ -7,13 +7,15 @@ version = "0.4.0"
 BioMarkovChains = "f861b655-cb5f-42ce-b66a-341b542d4f2c"
 BioSequences = "7e6ae17a-c86d-528c-b3b9-7f778a29fe59"
 FASTX = "c2308a5c-f048-11e8-3e8a-31650f418d12"
+GenomicFeatures = "899a7d2d-5c61-547b-bef9-6698a8d05446"
 IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 
 [compat]
 BioMarkovChains = "0.9"
 BioSequences = "3"
 FASTX = "2"
+GenomicFeatures = "3"
 IterTools = "1.4"
 PrecompileTools = "1"
 julia = "1"

src/GeneFinder.jl

@@ -25,6 +25,7 @@ using BioMarkovChains: BioMarkovChain, dnaseqprobability, ECOLICDS, ECOLINOCDS,
 using FASTX: FASTAReader, sequence
 using IterTools: takewhile, iterated
 using PrecompileTools: @setup_workload, @compile_workload
+using GenomicFeatures: GenomicFeatures, AbstractGenomicInterval, GenomicInterval, Strand, summary, groupname, leftposition, rightposition, strand, metadata, STRAND_POS
 
 include("algorithms/naivefinder.jl")
 include("types.jl")

src/algorithms/naivefinder.jl

@@ -1,4 +1,6 @@
-export naivefinder, naivefinderscored, log_odds_ratio_decision_rule, lordr
+export naivefinder, naivefinderscored, log_odds_ratio_decision_rule, lordr, lors
+
+const lors = log_odds_ratio_score
 
 """
     locationiterator(sequence::NucleicSeqOrView{DNAAlphabet{N}}; alternative_start::Bool=false) where {N}
@@ -48,20 +50,21 @@ function naivefinder(
     sequence::NucleicSeqOrView{DNAAlphabet{N}};
     alternative_start::Bool = false,
     min_len::Int64 = 6,
+    scheme::Union{Nothing, Function} = nothing,
     kwargs...
 ) where {N}
     seqlen = length(sequence)
     framedict = Dict(0 => 3, 1 => 1, 2 => 2)
     orfs = Vector{ORF}()
+    seqname = get_variable_name(sequence)
     for strand in ('+', '-')
         seq = strand == '-' ? reverse_complement(sequence) : sequence
-
         @inbounds for location in @views _locationiterator(seq; alternative_start)
             if length(location) >= min_len
                 frame = strand == '+' ? framedict[location.start % 3] : framedict[(seqlen - location.stop + 1) % 3]
                 start = strand == '+' ? location.start : seqlen - location.stop + 1
                 stop = start + length(location) - 1
-                push!(orfs, ORF(start:stop, strand, frame, 0.0))
+                push!(orfs, ORF(seqname, start:stop, strand, frame, scheme)) #NaiveFinder()
             end
         end
     end
@@ -91,12 +94,13 @@ function naivefinderscored(
     sequence::NucleicSeqOrView{DNAAlphabet{N}};
     alternative_start::Bool = false,
     min_len::Int64 = 6,
-    scoringscheme::BioMarkovChain = ECOLICDS, # ECOLINOCDS
+    model::BioMarkovChain = ECOLICDS,
     kwargs...
 ) where {N}
     seqlen = length(sequence)
     framedict = Dict(0 => 3, 1 => 1, 2 => 2)
     orfs = Vector{ORF}()
+    seqname = get_variable_name(sequence)
     for strand in ('+', '-')
         seq = strand == '-' ? reverse_complement(sequence) : sequence
 
@@ -106,8 +110,8 @@ function naivefinderscored(
                 start = strand == '+' ? location.start : seqlen - location.stop + 1
                 stop = start + length(location) - 1
                 # score = -10log10(dnaseqprobability(seq[start:stop], scoringscheme))
-                score = log_odds_ratio_score(seq[start:stop], scoringscheme)
-                push!(orfs, ORF(start:stop, strand, frame, score))
+                score = log_odds_ratio_score(seq[start:stop], model)
+                push!(orfs, ORF(seqname, start:stop, strand, frame, NaiveFinderScored(), lors, score))
             end
         end
     end
@@ -154,7 +158,7 @@ sequence = dna"ATGGCATCTAG"
 iscoding(sequence)  # Returns: true or false
 ```
 """
-function log_odds_ratio_decision_rule( #log_odds_ratio_decision lordr/cudr/kfdr/aadr
+function lordr( #log_odds_ratio_decision, also lordr/cudr/kfdr/aadr
     sequence::NucleicSeqOrView{DNAAlphabet{N}};
     codingmodel::BioMarkovChain = ECOLICDS,
     noncodingmodel::BioMarkovChain = ECOLINOCDS,
@@ -176,7 +180,7 @@ function log_odds_ratio_decision_rule( #log_odds_ratio_decision lordr/cudr/kfdr/
     end
 end
 
-const lordr = log_odds_ratio_decision_rule # criteria
+const log_odds_ratio_decision_rule = lordr # criteria
 
 ## Another alternative:
 

src/extended.jl

@@ -2,7 +2,9 @@
 import Base: isless, iterate, sort, getindex
 
 
-Base.isless(a::ORF, b::ORF) = isless(a.location, b.location)
+# Base.isless(a::ORF, b::ORF) = isless(a.location, b.location)
+Base.isless(a::ORF, b::ORF) = isless(a.first:a.last, b.first:b.last)
+# Base.isless(a::ORF, b::ORF) = isless(a.location, b.location)
 # Base.isless(a::ORF, b::ORF) = isless(a.score, b.score)
 
 # Base.sort(v::Vector{<:ORF}; kwargs...) = sort(v, by = _orf_sort_key)
@@ -12,11 +14,19 @@ Base.isless(a::ORF, b::ORF) = isless(a.location, b.location)
 #TODOs: how to make more robust the getindex method? confroning the frames?
 # Base.getindex(sequence::NucleicSeqOrView{A}, orf::ORF) where {A} = orf.strand == '+' ? (@view sequence[orf.location]) : reverse_complement(@view sequence[orf.location])
 
+# function getindex(sequence::NucleicSeqOrView{A}, orf::ORF) where {A}
+#     if orf.strand == '+'
+#         return @view sequence[orf.location]
+#     else
+#         return reverse_complement(@view sequence[orf.location])
+#     end
+# end
+
 function getindex(sequence::NucleicSeqOrView{A}, orf::ORF) where {A}
-    if orf.strand == '+'
-        return @view sequence[orf.location]
+    if orf.strand == '+' || orf.strand == STRAND_POS
+        return @view sequence[orf.first:orf.last]
     else
-        return reverse_complement(@view sequence[orf.location])
+        return reverse_complement(@view sequence[orf.first:orf.last])
     end
 end
 

src/findorfs.jl

@@ -54,6 +54,7 @@ function findorfs(
 end
 
 
+
 ### Some possible ideas:
 
 # function findorfs(