Merge pull request #16 from camilogarciabotero/dev

Exploring the log-odds ratio and other methods

CHANGELOG.md

@@ -6,6 +6,14 @@ 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)
 
+## [0.6.0]
+
+- Update docs version to `1.0`
+- Update docs logo
+- Add new methods for to calculate the log-odds ratio of a sequence and a model (#)
+- Extend new methods from the booleans/predicates of `DiscreteMarkovChains.jl` package
+- Clean code and make it more lightweight by using weakdeps
+
 ## [0.5.0]
 
 - Improve count of transitions 10X thanks to `VectorizedKmers.jl`

Project.toml

@@ -1,11 +1,17 @@
 name = "BioMarkovChains"
 uuid = "f861b655-cb5f-42ce-b66a-341b542d4f2c"
 authors = ["Camilo García<[email protected]>"]
-version = "0.5.0"
+version = "0.6.0"
+
+[weakdeps]
+DiscreteMarkovChains = "8abcb7ef-b365-4f7b-ac38-56893fb62f9f"
+MarkovChainHammer = "38c40fd0-bccb-4723-b30d-b2caea0ad8d9"
+
+[extensions]
+DiscreteMarkovChainsExt = "DiscreteMarkovChains"
 
 [deps]
 BioSequences = "7e6ae17a-c86d-528c-b3b9-7f778a29fe59"
-MarkovChainHammer = "38c40fd0-bccb-4723-b30d-b2caea0ad8d9"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 StatsAPI = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
@@ -15,14 +21,15 @@ VectorizedKmers = "2ef45bd6-4c2b-43d8-88b3-40597287359a"
 
 [compat]
 BioSequences = "3"
-MarkovChainHammer = "0.0.9, 0.0.10"
+DiscreteMarkovChains = "0.2"
+MarkovChainHammer = "0.0.11"
 PrecompileTools = "1"
 StatsAPI = "1.6"
 StatsBase = "0.34.0"
 TestItemRunner = "0.2"
 TestItems = "0.1"
-VectorizedKmers = "0.5, 0.6"
-julia = "1.9"
+VectorizedKmers = "0.7"
+julia = "1"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

docs/src/index.md

@@ -1,4 +1,10 @@
-![]("../assets/logo.svg")
+```@raw html
+
+<p align="center">
+  <img src="assets/logo.svg" height="150"><br/>
+  <i>Representing biological sequences as Markov chains</i>
+</p>
+```
 
 # BioMarkovChains
 

ext/DiscreteMarkovChainsExt.jl

@@ -0,0 +1,30 @@
+module DiscreteMarkovChainsExt
+
+using BioMarkovChains
+using DiscreteMarkovChains:
+    DiscreteMarkovChains, DiscreteMarkovChain,
+    is_absorbing, is_ergodic, is_regular, is_reversible
+
+## Boolean functions ##
+
+import DiscreteMarkovChains: is_absorbing, is_ergodic, is_regular, is_reversible
+
+function DiscreteMarkovChains.is_absorbing(bmc::BioMarkovChain)
+    return is_absorbing(DiscreteMarkovChain(bmc.tpm))
+end
+
+function DiscreteMarkovChains.is_ergodic(bmc::BioMarkovChain)
+    return is_ergodic(DiscreteMarkovChain(bmc.tpm))
+end
+
+function DiscreteMarkovChains.is_regular(bmc::BioMarkovChain)
+    return is_regular(DiscreteMarkovChain(bmc.tpm))
+end
+
+function DiscreteMarkovChains.is_reversible(bmc::BioMarkovChain)
+    return is_reversible(DiscreteMarkovChain(bmc.tpm))
+end
+
+## 
+
+end

ext/MarkovChainHammerExt.jl

@@ -0,0 +1,6 @@
+# module MarkovChainHammerExt
+
+# using BioMarkovChains
+# using MarkovChainHammer.Trajectory: generate
+
+# end

src/BioMarkovChains.jl

@@ -28,8 +28,6 @@ using BioSequences:
     alphabet,
     @biore_str
 
-
-using MarkovChainHammer.Trajectory: generate
 using PrecompileTools: @setup_workload, @compile_workload
 using TestItems: @testitem
 using StatsAPI: StatsAPI, fit, fit!
@@ -43,18 +41,18 @@ include("utils.jl")
 export transitions, hasprematurestop, randbmc
 
 include("transitions.jl")
-export initials, transition_count_matrix, transition_probability_matrix, dnaseqprobability
+export initials, transition_count_matrix, transition_probability_matrix, dnaseqprobability, logoddsratio, logoddsratioscore
 
 include("models.jl")
 export ECOLICDS, ECOLINOCDS
 
 include("perronfrobenius.jl")
 export perronfrobenius, generatedna
 
-include("extended.jl")
-
 include("constants.jl")
 
+include("extended.jl")
+
 @setup_workload begin
     # Putting some things in `@setup_workload` instead of `@compile_workload` can reduce the size of the
     # precompile file and potentially make loading faster.

src/extended.jl

@@ -1,4 +1,4 @@
-import Base: show, length
+import Base: show, length, eltype
 import StatsAPI: fit!
 
 function Base.show(io::IO, model::BioMarkovChain)
@@ -46,6 +46,10 @@ end
 
 Base.length(bmc::BioMarkovChain) = length(bmc.inits)
 
+function Base.eltype(bmc::BioMarkovChain)
+    return length(bmc) == 4 ? DNA : AminoAcid
+end
+
 """
     fit!(bmc::BMC, inits:Vector{Float64}, tpm::Matrix{Float64})
 
@@ -57,4 +61,4 @@ function StatsAPI.fit!(bmc::BMC, inits::Vector{Float64}, tpm::Matrix{Float64})
     bmc.tpm .= tpm
     foreach(sum_to_one!, eachrow(bmc.tpm))
     return nothing
-end
+end

src/models.jl

@@ -21,5 +21,33 @@ const ECOLINOCDS = begin
 
     inits = [0.262, 0.239, 0.240, 0.259]
 
+    BMC(tpm, inits)
+end
+
+
+const CPGPOS = begin
+    tpm = [
+        0.180 0.274 0.426 0.120
+        0.171 0.368 0.274 0.188
+        0.161 0.339 0.375 0.125 
+        0.079 0.355 0.384 0.182
+    ]
+
+    inits = [0.262, 0.239, 0.240, 0.259] # not stablished
+
+    BMC(tpm, inits)
+end
+
+
+const CPGNEG = begin
+    tpm = [
+        0.300 0.205 0.285 0.210 
+        0.322 0.298 0.078 0.302
+        0.248 0.246 0.298 0.208
+        0.177 0.239 0.292 0.292
+    ]
+
+    inits = [0.262, 0.239, 0.240, 0.259] # not stablished
+
     BMC(tpm, inits)
 end

src/perronfrobenius.jl

@@ -22,19 +22,15 @@ pf = perronfrobenius(sequence, n)
 """
 function perronfrobenius(sequence::LongNucOrView{4}; n::Int64=1)
     tpm = transition_probability_matrix(sequence, n)
-    pf = transpose(tpm)
-    return copy(pf)
+    return copy(tpm')
 end
 
 function perronfrobenius(bmc::BioMarkovChain)
-   tpm = bmc.tpm
-   pf = transpose(tpm)
-   return copy(pf)
+   return copy(bmc.tpm')
 end
 
 function perronfrobenius(tpm::Matrix{Float64})
-    pf = transpose(tpm)
-    return copy(pf)
+    return copy(tpm')
 end
 
 """
@@ -73,7 +69,7 @@ newdna = generatedna(pf, 100)
 function generatedna(pf::Matrix{Float64}, steps::Int64; extended_alphabet::Bool = false)
     newseq = LongDNA{4}()
     # pf = transpose(tpm) # The Perron-Frobenius matrix
-    trajectory = generate(pf, steps)
+    trajectory = MarkovChainHammer.generate(pf, steps)
     for i in trajectory
         newseq = append!(newseq, _int_to_dna(i; extended_alphabet))
     end