diff --git a/docs/doc.main b/docs/doc.main index eccd60f50df0..fa22649159fa 100644 --- a/docs/doc.main +++ b/docs/doc.main @@ -237,14 +237,15 @@ ], "Combinatorics" => [ - "Combinatorics/graphs.md", - "Combinatorics/matroids.md", - "Combinatorics/simplicialcomplexes.md", - "Enumerative combinatorics" => [ - "Combinatorics/EnumerativeCombinatorics/partitions.md", - "Combinatorics/EnumerativeCombinatorics/tableaux.md", - "Combinatorics/EnumerativeCombinatorics/schur_polynomials.md", - ], + "Combinatorics/graphs.md", + "Combinatorics/matroids.md", + "Combinatorics/simplicialcomplexes.md", + "Combinatorics/phylogenetic_trees.md", + "Enumerative combinatorics" => [ + "Combinatorics/EnumerativeCombinatorics/partitions.md", + "Combinatorics/EnumerativeCombinatorics/tableaux.md", + "Combinatorics/EnumerativeCombinatorics/schur_polynomials.md", + ], ], "Straight Line Programs" => [ diff --git a/docs/oscar_references.bib b/docs/oscar_references.bib index 29e7251a13fb..f14b6e563dd0 100644 --- a/docs/oscar_references.bib +++ b/docs/oscar_references.bib @@ -1721,6 +1721,16 @@ @Article{RSS03 year = {2003} } +@Book{SS03, + author = {Semple, Charles and Steel, Mike}, + title = {Phylogenetics}, + series = {Oxf. Lect. Ser. Math. Appl.}, + volume = {24}, + publisher = {Oxford University Press}, + year = {2003}, + fseries = {Oxford Lecture Series in Mathematics and its Applications} +} + @Article{SS12, author = {Savage, Carla D. and Schuster, Michael J.}, title = {Ehrhart series of lecture hall polytopes and Eulerian polynomials for inversion sequences}, diff --git a/docs/src/Combinatorics/phylogenetic_trees.md b/docs/src/Combinatorics/phylogenetic_trees.md new file mode 100644 index 000000000000..8058ddca29cd --- /dev/null +++ b/docs/src/Combinatorics/phylogenetic_trees.md @@ -0,0 +1,27 @@ +```@meta +CurrentModule = Oscar +``` + +# Phylogenetic Trees + +## Introduction + +Phylogenetic trees represent the evolutionary history of some species of consideration. +Here we consider phylogenetic trees with branch lengths as defined in [SS03](@cite). + +## Construction + +```@docs +phylogenetic_tree +``` + +## Some Helpful Functions + +```@docs +adjacency_tree +equidistant +cophenetic_matrix +taxa +newick +tropical_median_consensus +``` diff --git a/src/Combinatorics/Graphs/structs.jl b/src/Combinatorics/Graphs/structs.jl index d7964b862ced..e5b44e77610b 100644 --- a/src/Combinatorics/Graphs/structs.jl +++ b/src/Combinatorics/Graphs/structs.jl @@ -2,6 +2,6 @@ import Oscar: Polymake import Oscar.Polymake: Directed, Undirected struct Graph{T <: Union{Directed, Undirected}} - pm_graph::Polymake.Graph{T} + pm_graph::Polymake.Graph{T} end diff --git a/src/Combinatorics/PhylogeneticTrees.jl b/src/Combinatorics/PhylogeneticTrees.jl new file mode 100644 index 000000000000..5a5a17a9e08a --- /dev/null +++ b/src/Combinatorics/PhylogeneticTrees.jl @@ -0,0 +1,284 @@ +struct PhylogeneticTree{T <: Union{Float64, QQFieldElem}} + pm_ptree::Polymake.LibPolymake.BigObjectAllocated +end + +function pm_object(PT::PhylogeneticTree) + return PT.pm_ptree +end + +@doc raw""" + phylogenetic_tree(T::Type{<:Union{Float64, QQFieldElem}}, newick::String) + +Constructs a phylogenetic tree with Newick representation `newick`. `T` indicates +the numerical type of the edge lengths. + +# Examples +Make a phylogenetic tree with 4 leaves from its Newick representation and print +its taxa and cophenetic matrix. +```jldoctest +julia> phylo_t = phylogenetic_tree(Float64, "((H:3,(C:1,B:1):2):1,G:4);"); + +julia> taxa(phylo_t) +4-element Vector{String}: + "B" + "C" + "G" + "H" + +julia> cophenetic_matrix(phylo_t) +4×4 Matrix{Float64}: + 0.0 2.0 8.0 6.0 + 2.0 0.0 8.0 6.0 + 8.0 8.0 0.0 8.0 + 6.0 6.0 8.0 0.0 +``` +""" +function phylogenetic_tree(T::Type{<:Union{Float64, QQFieldElem}}, newick::String) + pm_ptree = Polymake.graph.PhylogeneticTree{Polymake.convert_to_pm_type(T)}(NEWICK = newick) + + # load graph properties + pm_ptree.ADJACENCY + + return PhylogeneticTree{T}(pm_ptree) +end + +@doc raw""" + phylogenetic_tree(M::Matrix{T}, taxa::Vector{String}) where T <: Union{Float64, QQFieldElem} + +Constructs a phylogenetic tree with cophenetic matrix `M` and taxa `taxa`. The matrix `M` must be +ultrametric, otherwise an error will be thrown. + +# Examples +Make a phylogenetic tree on 4 taxa with given cophenetic matrix and print one Newick representation. + +```jldoctest +julia> mat = [0. 2 8 6; 2 0 8 6; 8 8 0 8; 6 6 8 0] +4×4 Matrix{Float64}: + 0.0 2.0 8.0 6.0 + 2.0 0.0 8.0 6.0 + 8.0 8.0 0.0 8.0 + 6.0 6.0 8.0 0.0 + +julia> tax = ["Bonobo", "Chimpanzee", "Gorilla", "Human"] +4-element Vector{String}: + "Bonobo" + "Chimpanzee" + "Gorilla" + "Human" + +julia> tree_mat = phylogenetic_tree(mat, tax); + +julia> newick(tree_mat) +"Gorilla:4,(Human:3,(Bonobo:1,Chimpanzee:1):2):1;" +``` +""" +function phylogenetic_tree(M::Matrix{Float64}, taxa::Vector{String}) + n_taxa = length(taxa) + @req (n_taxa, n_taxa) == size(M) "Number of taxa should match the rows and columns of the given matrix" + pm_ptree = Polymake.graph.PhylogeneticTree{Float64}(COPHENETIC_MATRIX = M, TAXA = taxa) + return PhylogeneticTree{Float64}(pm_ptree) +end + +function phylogenetic_tree(M::QQMatrix, taxa::Vector{String}) + n_taxa = length(taxa) + @req (n_taxa, n_taxa) == size(M) "Number of taxa should match the rows and columns of the given matrix" + pm_ptree = Polymake.graph.PhylogeneticTree{Rational}( + COPHENETIC_MATRIX = M, TAXA = taxa + ) + return PhylogeneticTree{QQFieldElem}(pm_ptree) +end + +@doc raw""" + adjacency_tree(ptree::PhylogeneticTree) + +Returns the underlying graph of the phylogenetic tree `ptree`. + +# Examples +Make a phylogenetic tree with given Newick format and print its underlying graph. + +```jldoctest +julia> ptree = phylogenetic_tree(Float64, "((H:3,(C:1,B:1):2):1,G:4);"); + +julia> adjacency_tree(ptree) +Undirected graph with 7 nodes and the following edges: +(2, 1)(3, 2)(4, 2)(5, 4)(6, 4)(7, 1) +``` +""" +function adjacency_tree(ptree::PhylogeneticTree) + return Graph{Undirected}(pm_object(ptree).ADJACENCY) +end + +function Base.show(io::IO, ptree::PhylogeneticTree{T}) where T + print(io, "Phylogenetic tree with $T type coefficients") +end + +@doc raw""" + equidistant(ptree::PhylogeneticTree) + +Checks if the phylogenetic tree `ptree` is equidistant. + +# Examples +Make a phylogenetic tree with given Newick format and check if it is equidistant. + +```jldoctest +julia> ptree = phylogenetic_tree(Float64, "((H:3,(C:1,B:1):2):1,G:4);"); + +julia> equidistant(ptree) +true +``` +""" +function equidistant(ptree::PhylogeneticTree) + return pm_object(ptree).EQUIDISTANT::Bool +end + + +@doc raw""" + cophenetic_matrix(ptree::PhylogeneticTree) + +Returns the cophenetic matrix of the phylogenetic tree `ptree`. + +# Examples +Make a phylogenetic tree with given Newick format and print its cophenetic matrix. + +```jldoctest +julia> ptree = phylogenetic_tree(Float64, "((H:3,(C:1,B:1):2):1,G:4);"); + +julia> cophenetic_matrix(ptree) +4×4 Matrix{Float64}: + 0.0 2.0 8.0 6.0 + 2.0 0.0 8.0 6.0 + 8.0 8.0 0.0 8.0 + 6.0 6.0 8.0 0.0 +``` +""" +function cophenetic_matrix(ptree::PhylogeneticTree{Float64}) + return convert(Matrix, pm_object(ptree).COPHENETIC_MATRIX)::Matrix{Float64} +end + +function cophenetic_matrix(ptree::PhylogeneticTree{QQFieldElem}) + return matrix(QQ, pm_object(ptree).COPHENETIC_MATRIX)::QQMatrix +end + +@doc raw""" + taxa(ptree::PhylogeneticTree) + +Returns the taxa of the phylogenetic tree `ptree`. + +# Examples +Make a phylogenetic tree with given Newick format and print its taxa. + +```jldoctest +julia> ptree = phylogenetic_tree(Float64, "((H:3,(C:1,B:1):2):1,G:4);"); + +julia> taxa(ptree) +4-element Vector{String}: + "B" + "C" + "G" + "H" +``` +""" +function taxa(ptree::PhylogeneticTree) + return convert(Array{String}, pm_object(ptree).TAXA)::Array{String} +end + +@doc raw""" + newick(ptree::PhylogeneticTree) + +Returns a Newick representation of the phylogenetic tree `ptree`. + +# Examples +Make a phylogenetic tree from a matrix and print a Newick representation of it. + +```jldoctest +julia> mat = [0. 2 8 6; 2 0 8 6; 8 8 0 8; 6 6 8 0] +4×4 Matrix{Float64}: + 0.0 2.0 8.0 6.0 + 2.0 0.0 8.0 6.0 + 8.0 8.0 0.0 8.0 + 6.0 6.0 8.0 0.0 + +julia> tax = ["Bonobo", "Chimpanzee", "Gorilla", "Human"] +4-element Vector{String}: + "Bonobo" + "Chimpanzee" + "Gorilla" + "Human" + +julia> tree_mat = phylogenetic_tree(mat, tax); + +julia> newick(tree_mat) +"Gorilla:4,(Human:3,(Bonobo:1,Chimpanzee:1):2):1;" +``` +""" +function newick(ptree::PhylogeneticTree) + return convert(String, pm_object(ptree).NEWICK)::String +end + + +@doc raw""" + tropical_median_consensus(arr::Vector{PhylogeneticTree{T}}) + +Computes the tropical median consensus tree of the phylogenetic trees from +the vector `arr`. + +# Examples +Compute the tropical median consensus of three trees and print one of its +Newick representations. + +```jldoctest +julia> t1 = phylogenetic_tree(Float64, "((H:30,(C:10,B:10):20):10,G:40);"); + +julia> t2 = phylogenetic_tree(Float64, "(((H:10,C:10):20,B:30):10,G:40);"); + +julia> t3 = phylogenetic_tree(Float64, "((H:25,C:25):15,(B:15,G:15):25);"); + +julia> arr = [t1, t2, t3]; + +julia> tc = tropical_median_consensus(arr); + +julia> newick(tc) +"G:40,(B:35,(C:30,H:30):5):5;" +``` +""" +function tropical_median_consensus(arr::Vector{PhylogeneticTree{T}}) where {T <: Union{Float64, QQFieldElem}} + + n = length(arr) + @req n > 0 "The vector must not be empty" + + phylo_type = Polymake.bigobject_type(pm_object(first(arr))) + pm_arr = Polymake.Array{Polymake.BigObject}(phylo_type, n) + + pm_arr .= pm_object.(arr) + + pm_cons_tree = Polymake.tropical.tropical_median_consensus(pm_arr) + return PhylogeneticTree{T}(pm_cons_tree) +end + + +@doc raw""" + tropical_median_consensus(trees::Vararg{PhylogeneticTree, N}) where {N} + +Computes the tropical median consensus tree of any number of phylogenetic trees +given as parameters. + +# Examples +Compute the tropical median consensus of three trees and print one of its +Newick representations. + +```jldoctest +julia> t1 = phylogenetic_tree(Float64, "((H:30,(C:10,B:10):20):10,G:40);"); + +julia> t2 = phylogenetic_tree(Float64, "(((H:10,C:10):20,B:30):10,G:40);"); + +julia> t3 = phylogenetic_tree(Float64, "((H:25,C:25):15,(B:15,G:15):25);"); + +julia> tc = tropical_median_consensus(t1, t2, t3); + +julia> newick(tc) +"G:40,(B:35,(C:30,H:30):5):5;" +``` +""" +function tropical_median_consensus(trees::Vararg{PhylogeneticTree, N}) where {N} + return tropical_median_consensus(collect(trees)) +end diff --git a/src/Oscar.jl b/src/Oscar.jl index 2c02dfa3a82d..f91d5e4deded 100644 --- a/src/Oscar.jl +++ b/src/Oscar.jl @@ -239,6 +239,7 @@ include("Combinatorics/OrderedMultiIndex.jl") include("Combinatorics/Matroids/JMatroids.jl") include("Combinatorics/Compositions.jl") include("Combinatorics/EnumerativeCombinatorics/EnumerativeCombinatorics.jl") +include("Combinatorics/PhylogeneticTrees.jl") include("StraightLinePrograms/StraightLinePrograms.jl") include("Rings/lazypolys.jl") # uses StraightLinePrograms diff --git a/src/Serialization/Combinatorics.jl b/src/Serialization/Combinatorics.jl index 67f49b2f7ac6..aea58af8e125 100644 --- a/src/Serialization/Combinatorics.jl +++ b/src/Serialization/Combinatorics.jl @@ -49,3 +49,17 @@ function load_object(s::DeserializerState, K::Type{SimplicialComplex}) bigobject = Polymake.call_function(:common, :deserialize_json_string, json(s.obj)) return K(bigobject) end + +############################################################################### +## Phylogenetic Trees +############################################################################### +@register_serialization_type PhylogeneticTree + +function save_object(s::SerializerState, PT::PhylogeneticTree) + save_object(s, pm_object(PT)) +end + +function load_object(s::DeserializerState, T::Type{PhylogeneticTree}, dict::Dict) + bigobject = Polymake.call_function(:common, :deserialize_json_string, json(dict)) + return T(bigobject) +end diff --git a/src/Serialization/polymake.jl b/src/Serialization/polymake.jl index d6de313dde0e..59246167da25 100644 --- a/src/Serialization/polymake.jl +++ b/src/Serialization/polymake.jl @@ -12,6 +12,10 @@ function load_from_polymake(::Type{Graph{T}}, jsondict::Dict{Symbol, Any}) where return Graph{T}(polymake_object) end +function load_from_polymake(::Type{PhylogeneticTree{T}}, jsondict::Dict{Symbol, Any}) where {T <: Union{Float64, Int, QQFieldElem}} + polymake_object = Polymake.call_function(:common, :deserialize_json_string, json(jsondict)) + return PhylogeneticTree{T}(polymake_object) +end const polymake2OscarTypes = Dict{String, Type}([ "polytope::Cone" => Cone{QQFieldElem}, @@ -22,6 +26,8 @@ const polymake2OscarTypes = Dict{String, Type}([ "topaz::SimplicialComplex" => SimplicialComplex, "common::GraphAdjacency" => Graph{Undirected}, "common::GraphAdjacency" => Graph{Directed}, + "graph::PhylogeneticTree" => PhylogeneticTree{QQFieldElem}, + "graph::PhylogeneticTree" => PhylogeneticTree{Float64} ]) @register_serialization_type Polymake.BigObjectAllocated "Polymake.BigObject" uses_id diff --git a/src/exports.jl b/src/exports.jl index 1678978308fa..a4533346a5b6 100644 --- a/src/exports.jl +++ b/src/exports.jl @@ -188,6 +188,7 @@ export add_edge! export add_gluing! export add_vertex! export add_vertices! +export adjacency_tree export adjacent_chamber export adjoint_ideal export affine_algebra @@ -389,6 +390,7 @@ export coordinate_names_of_torus export coordinate_ring export coordinate_ring_of_torus export coordinates +export cophenetic_matrix export corank export core export corresponding_bilinear_form @@ -494,6 +496,7 @@ export equidimensional_decomposition_radical export equidimensional_decomposition_weak export equidimensional_hull export equidimensional_hull_radical +export equidistant export euler_characteristic export euler_phi export expand @@ -1043,6 +1046,7 @@ export neglex export negwdeglex export negwdegrevlex export neighbors +export newick export newton_polytope export nilpotency_class, has_nilpotency_class, set_nilpotency_class export noether_normalization @@ -1145,6 +1149,7 @@ export permutation_group export permutation_matrix export permutation_of_terms export permuted +export phylogenetic_tree export picard_class export picard_group export picard_index @@ -1410,6 +1415,7 @@ export syzygy_generators export tail export tangent_lines export tangent_space +export taxa export tensor_product export terms export tetrahedron @@ -1443,6 +1449,7 @@ export trivial_divisor_class export trivial_morphism export trivial_subgroup, has_trivial_subgroup, set_trivial_subgroup export tropical_matrix +export tropical_median_consensus export tropical_pluecker_vector export tropical_polynomial export tropical_variety diff --git a/test/Combinatorics/PhylogeneticTrees.jl b/test/Combinatorics/PhylogeneticTrees.jl new file mode 100644 index 000000000000..e0221dc814ea --- /dev/null +++ b/test/Combinatorics/PhylogeneticTrees.jl @@ -0,0 +1,19 @@ +@testset "Phylogenetic Trees" begin + ptree1 = phylogenetic_tree(Float64, "((A:4,B:4):5,C:9);") + test_tree = graph_from_edges([[2, 1], [3, 2], [4, 2], [5, 1]]) + @test is_isomorphic(adjacency_tree(ptree1), test_tree) + @test equidistant(ptree1) == true + + ptree1_test = phylogenetic_tree([0 4 9; 4 0 9.0; 9 9 0], ["A", "B", "C"]) + tmc = tropical_median_consensus(ptree1_test, ptree1) + tmc_v = tropical_median_consensus([ptree1_test, ptree1]) + @test newick(tmc) == "C:9,(A:6.5,B:6.5):2.5;" + @test newick(tmc_v) == "C:9,(A:6.5,B:6.5):2.5;" + + M = matrix(QQ, [0 4 9; 4 0 9; 9 9 0]) + ptree2 = phylogenetic_tree(M, ["a", "b", "c"]) + + @test cophenetic_matrix(ptree2) == matrix(QQ, [0 4 9; 4 0 9; 9 9 0]) + @test newick(ptree2) == "c:9/2,(a:2,b:2):5/2;" + @test taxa(ptree2) == ["a", "b", "c"] +end