Phylogenetics
The Phylogenetics module in OSCAR provides a comprehensive toolkit for the algebraic and geometric study of phylogenetic models, specifically nucleotide substitution Markov processes on phylogenetic trees. It allows for the construction of various evolutionary models on phylogenetic trees, the creation of their associated polynomial rings and parametrizations, and the computation of their algebraic properties, such as phylogenetic invariants (polynomials in the venishing ideal of the model).
Model Constructors
Two main types of models can be constructed: a general PhylogeneticModel and a more specialized GroupBasedPhylogeneticModel.
PhylogeneticModel
A PhylogeneticModel is defined by a directed tree, a base field, the symbolic structure of the transition matrices for each edge, and a probability distribution at the root.
PhylogeneticModel — Type
PhylogeneticModel{GT, L, M, R}A data structure representing a general phylogenetic model on a directed tree.
This model defines a probability distribution on the states of the leaves of a phylogenetic tree, based on a root distribution and transition matrices on the edges. The model is defined over a base field (e.g., QQ for rational numbers).
Type Parameters:
GT: The graph type, expected to beGraph{Directed},PhylogeneticTreeorPhylogeneticNetwork.L: The type of the graph labelings, aNamedTuple.M: The type of elements in the transition matrix structure, typically aVarName(aSymbol) or a polynomial ring element (MPolyRingElem).R: The type of elements in the root distribution vector.
Fields:
base_field::Field: The base field of the model's parameters.graph::GT: The underlying directed graph (tree) of the model.labelings::L: A named tuple of graph maps for internal use.trans_matrix_structure::Matrix{M}: A matrix representing the symbolic structure of the transition matrices for each edge.root_distribution::Vector{T}: A vector representing the distribution at the root of the tree.n_states::Int: The number of states in the model.model_parameter_name::VarName: The symbolic name for the variables of the model ring.
This function is part of the experimental code in Oscar. Please read here for more details.
phylogenetic_model — Function
phylogenetic_model(F::Field, G::AbstractGraph{Directed}, trans_matrix_structure::Matrix, root_distribution::Union{Nothing, Vector} = nothing, varnames::VarName="p")
phylogenetic_model(G::AbstractGraph{Directed}, trans_matrix_structure::Matrix, root_distribution::Union{Nothing, Vector} = nothing, varnames::VarName="p")
phylogenetic_model(G::AbstractGraph{Directed}, trans_matrix_structure::Matrix{M}, root_distribution::Vector{T}, varnames::VarName="p") where {M <: MPolyRing{<:MPolyRingElem}, R <: MPolyRing}Construct a PhylogeneticModel from a field F, a graph G, a symbolic transition matrix trans_matrix_structure, and an optional root_distribution.
G can either be of type Graph{Directed}, PhylogeneticTree or PhylogeneticNetwork.
If root_distribution is not provided, a uniform distribution is assumed. If F is not provided, it constructs a PhylogeneticModel using the default rational field (QQ).
If trans_matrix_structure is a Matrix{<: MPolyRing{<:MPolyRingElem}} and root_distribution is a Vector{<: MPolyRing}, then the constructor ensures the rings of the root distribution parameters and the transition matrices are compatible.
Examples
julia> tree = phylogenetic_tree(graph_from_edges(Directed,[[4,1], [4,2], [4,3]]));
julia> M = [:m11 :m12; :m21 :m22];
julia> root_dist = [:r1, :r2];
julia> PM = phylogenetic_model(tree, M, root_dist)
Phylogenetic model on a tree with 3 leaves and 3 edges
with root distribution [r1, r2] and transition matrices of the form
[:m11 :m12;
:m21 :m22].This function is part of the experimental code in Oscar. Please read here for more details.
phylogenetic_model(PM::GroupBasedPhylogeneticModel)Return the underlying PhylogeneticModel.
Examples
julia> tree = phylogenetic_tree(graph_from_edges(Directed,[[4,1], [4,2], [4,3]]));
julia> PM = jukes_cantor_model(tree);
julia> phylogenetic_model(PM)
Phylogenetic model on a tree with 3 leaves and 3 edges
with root distribution [1//4, 1//4, 1//4, 1//4] and transition matrices of the form
[:a :b :b :b;
:b :a :b :b;
:b :b :a :b;
:b :b :b :a]. This function is part of the experimental code in Oscar. Please read here for more details.
Here are some examples of how to construct phylogenetic models. The Jukes Cantor model on a 3-leaf tree can be constructed:
julia> tree = graph_from_edges(Directed,[[4,1], [4,2], [4,3]]);
julia> M_JC = [:a :b :b :b;
:b :a :b :b;
:b :b :a :b;
:b :b :b :a];
julia> PM_jukesCantor = phylogenetic_model(tree, M_JC)
Phylogenetic model on a tree with 3 leaves and 3 edges
with root distribution [1//4, 1//4, 1//4, 1//4] and transition matrices of the form
[:a :b :b :b;
:b :a :b :b;
:b :b :a :b;
:b :b :b :a].And the Tamura-Nei model (TN93) in the same tree can be constructed as follows:
julia> tree = graph_from_edges(Directed,[[4,1], [4,2], [4,3]]);
julia> Rp, r = rational_function_field(QQ, :r => 1:4);
julia> RM, (a, b, c, d) = polynomial_ring(Rp, [:a, :b, :c, :d]);
julia> M_TN93 = [a*r[1] c*r[2] b*r[3] b*r[4];
c*r[1] a*r[2] b*r[3] b*r[4];
b*r[1] b*r[2] a*r[3] d*r[4];
b*r[1] b*r[2] d*r[3] a*r[4]];
julia> PM_TN93 = PhylogeneticModel(tree, M_TN93, r)
Phylogenetic model on a tree with 3 leaves and 3 edges
with root distribution [r[1], r[2], r[3], r[4]] and transition matrices of the form
[r[1]*a r[2]*c r[3]*b r[4]*b;
r[1]*c r[2]*a r[3]*b r[4]*b;
r[1]*b r[2]*b r[3]*a r[4]*d;
r[1]*b r[2]*b r[3]*d r[4]*a]. GroupBasedPhylogeneticModel
For models exhibiting symmetries that can be captured by a finite abelian group, the GroupBasedPhylogeneticModel is used. This structure requires the symbolic Fourier parameters and the associated group, in addition to the standard phylogenetic model components.
GroupBasedPhylogeneticModel — Type
GroupBasedPhylogeneticModel{GT, L} <: GraphicalModel{GT, L}A data structure representing a group-based phylogenetic model.
This model is a specialization of PhylogeneticModel where the transition matrices and leaf probabilities can be described by a set of Fourier parameters related to a group structure.
Type Parameters:
GT: The graph type, expected to beGraph{Directed}.L: The type of the graph labelings.
Fields:
phylo_model::PhylogeneticModel: The underlying standard phylogenetic model.fourier_param_structure::Vector{<: VarName}: A vector representing the symbolic structure of the Fourier parameters.group::Vector{FinGenAbGroupElem}: The elements of the finite abelian group used for the model.model_parameter_name::VarName: The symbolic name for the model's Fourier parameters.
This function is part of the experimental code in Oscar. Please read here for more details.
group_based_phylogenetic_model — Function
group_based_phylogenetic_model(F::Field, G::AbstractGraph{Directed}, trans_matrix_structure::Matrix{<: VarName}, fourier_param_structure::Vector{<: VarName}, group::Union{Nothing, Vector{FinGenAbGroupElem}} = nothing, root_distribution::Union{Nothing, Vector} = nothing, varnames_phylo_model::VarName="p", varnames_group_based::VarName="q"))
group_based_phylogenetic_model(G::AbstractGraph{Directed}, trans_matrix_structure::Matrix{<: VarName}, fourier_param_structure::Vector{<: VarName}, group::Union{Nothing, Vector{FinGenAbGroupElem}} = nothing, root_distribution::Union{Nothing, Vector} = nothing, varnames_phylo_model::VarName="p", varnames_group_based::VarName="q")Construct a GroupBasedPhylogeneticModel from a field F, a graph G (of type Graph{Directed}, PhylogeneticTree or PhylogeneticNetwork), a symbolic transition matrix, symbolic Fourier parameters, an optional group, and optional root distribution.
If the group is not provided, a default is used based on the number of states: the group $Z_2$ for 2-state models, and the Klein four-group ($Z_2 \times Z_2$) for 4-state models. For other numbers of states, a group must be specified. If the root_distribution is not provided, a uniform distribution is assumed. If F is not provided, the model is constructed over the default rational field (QQ).
Examples
julia> tree = phylogenetic_tree(graph_from_edges(Directed,[[4,1], [4,2], [4,3]]));
julia> M = [:a :b; :b :a];
julia> fourier_params = [:x, :y];
julia> PM = group_based_phylogenetic_model(tree, M, fourier_params)
Group-based phylogenetic model on a tree with 3 leaves and 3 edges
with root distribution [1//2, 1//2],
transition matrices of the form
[:a :b;
:b :a]
and fourier parameters of the form [:x, :y].This function is part of the experimental code in Oscar. Please read here for more details.
For example, the Jukes-Cantor model can be defined as a group-based model:
julia> tree = graph_from_edges(Directed,[[4,1], [4,2], [4,3]]);
julia> M = [:a :b :b :b;
:b :a :b :b;
:b :b :a :b;
:b :b :b :a];
julia> x = [:x, :y, :y, :y];
julia> GroupBasedPhylogeneticModel(tree, M, x)
Group-based phylogenetic model on a tree with 3 leaves and 3 edges
with root distribution [1//4, 1//4, 1//4, 1//4],
transition matrices of the form
[:a :b :b :b;
:b :a :b :b;
:b :b :a :b;
:b :b :b :a]
and fourier parameters of the form [:x, :y, :y, :y].Classic phylogenetic models
Several classic phylogenetic models are pre-defined for convenience.
cavender_farris_neyman_model — Function
cavender_farris_neyman_model(F::Field, G::AbstractGraph{Directed})
cavender_farris_neyman_model(G::AbstractGraph{Directed})Constructs a GroupBasedPhylogeneticModel corresponding to the Cavender-Farris-Neyman model, a 2-state model with a single transition probability per edge. If the field F is not provided, the model is constructed over the default rational field (QQ).
Examples
julia> cavender_farris_neyman_model(graph_from_edges(Directed,[[4,1], [4,2], [4,3]]))
Group-based phylogenetic model on a tree with 3 leaves and 3 edges
with root distribution [1//2, 1//2],
transition matrices of the form
[:a :b;
:b :a]
and fourier parameters of the form [:x, :y].This function is part of the experimental code in Oscar. Please read here for more details.
jukes_cantor_model — Function
jukes_cantor_model(F::Field, G::AbstractGraph{Directed})
jukes_cantor_model(G::AbstractGraph{Directed})Constructs a GroupBasedPhylogeneticModel corresponding to the Jukes-Cantor model, a 4-state model where all non-diagonal transition probabilities are equal. If the field F is not provided, the model is constructed over the default rational field (QQ).
Examples
julia> tree = phylogenetic_tree(graph_from_edges(Directed,[[4,1], [4,2], [4,3]]));
julia> jukes_cantor_model(tree)
Group-based phylogenetic model on a tree with 3 leaves and 3 edges
with root distribution [1//4, 1//4, 1//4, 1//4],
transition matrices of the form
[:a :b :b :b;
:b :a :b :b;
:b :b :a :b;
:b :b :b :a]
and fourier parameters of the form [:x, :y, :y, :y].This function is part of the experimental code in Oscar. Please read here for more details.
kimura2_model — Function
kimura2_model(F::Field, G::AbstractGraph{Directed})
kimura2_model(G::AbstractGraph{Directed})Constructs a GroupBasedPhylogeneticModel corresponding to the Kimura 2-parameter model, a 4-state model with two non-diagonal transition probabilities. If the field F is not provided, the model is constructed over the default rational field (QQ).
Examples
julia> tree = phylogenetic_tree(graph_from_edges(Directed,[[4,1], [4,2], [4,3]]));
julia> kimura2_model(tree)
Group-based phylogenetic model on a tree with 3 leaves and 3 edges
with root distribution [1//4, 1//4, 1//4, 1//4],
transition matrices of the form
[:a :b :c :b;
:b :a :b :c;
:c :b :a :b;
:b :c :b :a]
and fourier parameters of the form [:x, :y, :z, :z].This function is part of the experimental code in Oscar. Please read here for more details.
kimura3_model — Function
kimura3_model(F::Field, G::AbstractGraph{Directed})
kimura3_model(G::AbstractGraph{Directed})Constructs a GroupBasedPhylogeneticModel corresponding to the Kimura 3-parameter model, a 4-state model with three non-diagonal different transition probabilities. If the field F is not provided, the model is constructed over the default rational field (QQ).
Examples
julia> tree = phylogenetic_tree(graph_from_edges(Directed,[[4,1], [4,2], [4,3]]));
julia> kimura3_model(tree)
Group-based phylogenetic model on a tree with 3 leaves and 3 edges
with root distribution [1//4, 1//4, 1//4, 1//4],
transition matrices of the form
[:a :b :c :d;
:b :a :d :c;
:c :d :a :b;
:d :c :b :a]
and fourier parameters of the form [:x, :y, :z, :t].This function is part of the experimental code in Oscar. Please read here for more details.
general_markov_model — Function
general_markov_model(F::Field, G::AbstractGraph{Directed})
general_markov_model(G::AbstractGraph{Directed})Constructs a PhylogeneticModel corresponding to a general Markov model with four states. All transition matrix entries and root distribution entries are treated as independent parameters. If the field F is not provided, the model is constructed over the default rational field (QQ).
Examples
julia> tree = phylogenetic_tree(graph_from_edges(Directed,[[4,1], [4,2], [4,3]]));
julia> general_markov_model(tree)
Phylogenetic model on a tree with 3 leaves and 3 edges
with root distribution [π1, π2, π3, π4] and transition matrices of the form
[:m11 :m12 :m13 :m14;
:m21 :m22 :m23 :m24;
:m31 :m32 :m33 :m34;
:m41 :m42 :m43 :m44].This function is part of the experimental code in Oscar. Please read here for more details.
general_time_reversible_model — Function
general_time_reversible_model(F::Field, G::AbstractGraph{Directed})
general_time_reversible_model(G::AbstractGraph{Directed})Constructs a PhylogeneticModel corresponding to a general time reversible model with four states. All transition matrix entries and root distribution entries are treated as independent parameters. If the field F is not provided, the model is constructed over the default rational field (QQ).
Examples
julia> tree = phylogenetic_tree(graph_from_edges(Directed,[[4,1], [4,2], [4,3]]));
julia> general_time_reversible_model(tree)
Phylogenetic model on a tree with 3 leaves and 3 edges
with root distribution [r[1], r[2], r[3], r[4]] and transition matrices of the form
[r[1]*a r[2]*b r[3]*c r[4]*d;
r[1]*b r[2]*a r[3]*e r[4]*f;
r[1]*c r[2]*e r[3]*a r[4]*g;
r[1]*d r[2]*f r[3]*g r[4]*a].This function is part of the experimental code in Oscar. Please read here for more details.
Properties
Objects of type GroupBasedPhylogeneticModel are built upon a PhylogeneticModel, which can be accessed directly:
phylogenetic_model — Method
phylogenetic_model(PM::GroupBasedPhylogeneticModel)Return the underlying PhylogeneticModel.
Examples
julia> tree = phylogenetic_tree(graph_from_edges(Directed,[[4,1], [4,2], [4,3]]));
julia> PM = jukes_cantor_model(tree);
julia> phylogenetic_model(PM)
Phylogenetic model on a tree with 3 leaves and 3 edges
with root distribution [1//4, 1//4, 1//4, 1//4] and transition matrices of the form
[:a :b :b :b;
:b :a :b :b;
:b :b :a :b;
:b :b :b :a]. This function is part of the experimental code in Oscar. Please read here for more details.
Other properties of a PhylogeneticModel and a GroupBasedPhylogeneticModel can be accessed using the following functions.
base_field — Method
base_field(PM::PhylogeneticModel)
base_field(PM::GroupBasedPhylogeneticModel)Return the base field of the underlying phylogenetic model.
Example
julia> tree = phylogenetic_tree(graph_from_edges(Directed,[[4,1], [4,2], [4,3]]));
julia> base_field(general_markov_model(tree))
Rational field
julia> base_field(jukes_cantor_model(tree))
Rational fieldThis function is part of the experimental code in Oscar. Please read here for more details.
n_states — Function
n_states(PM::PhylogeneticModel)
n_states(PM::GroupBasedPhylogeneticModel)Return the number of states of the underlying phylogenetic model.
Examples
julia> tree = phylogenetic_tree(graph_from_edges(Directed,[[4,1], [4,2], [4,3]]));
julia> n_states(general_markov_model(tree))
4
julia> n_states(cavender_farris_neyman_model(tree))
2This function is part of the experimental code in Oscar. Please read here for more details.
root_distribution — Function
root_distribution(PM::PhylogeneticModel)
root_distribution(PM::GroupBasedPhylogeneticModel)Return the root distribution vector of the underlying phylogenetic model.
Examples
julia> tree = phylogenetic_tree(graph_from_edges(Directed,[[4,1], [4,2], [4,3]]));
julia> root_distribution(general_markov_model(tree))
4-element Vector{Symbol}:
:π1
:π2
:π3
:π4
julia> root_distribution(cavender_farris_neyman_model(tree))
2-element Vector{QQFieldElem}:
1//2
1//2This function is part of the experimental code in Oscar. Please read here for more details.
transition_matrix — Function
transition_matrix(PM::PhylogeneticModel)
transition_matrix(PM::GroupBasedPhylogeneticModel)Return the structure of the transition matrices of the underlying phylogenetic model.
Examples
julia> tree = phylogenetic_tree(graph_from_edges(Directed,[[4,1], [4,2], [4,3]]));
julia> transition_matrix(general_markov_model(tree))
4×4 Matrix{Symbol}:
:m11 :m12 :m13 :m14
:m21 :m22 :m23 :m24
:m31 :m32 :m33 :m34
:m41 :m42 :m43 :m44
julia> transition_matrix(cavender_farris_neyman_model(tree))
2×2 Matrix{Symbol}:
:a :b
:b :aThis function is part of the experimental code in Oscar. Please read here for more details.
fourier_parameters — Function
fourier_parameters(PM::GroupBasedPhylogeneticModel)Return the symbolic structure of the Fourier parameters.
Examples
julia> tree = phylogenetic_tree(graph_from_edges(Directed,[[4,1], [4,2], [4,3]]));
julia> PM = jukes_cantor_model(tree);
julia> fourier_parameters(PM)
4-element Vector{Symbol}:
:x
:y
:y
:yThis function is part of the experimental code in Oscar. Please read here for more details.
group — Method
group(PM::GroupBasedPhylogeneticModel)Return the finite abelian group associated with the model.
Examples
julia> tree = phylogenetic_tree(graph_from_edges(Directed,[[4,1], [4,2], [4,3]]));
julia> PM = jukes_cantor_model(tree);
julia> group(PM)
4-element Vector{FinGenAbGroupElem}:
[0, 0]
[0, 1]
[1, 0]
[1, 1]This function is part of the experimental code in Oscar. Please read here for more details.
varnames — Method
varnames(PM::PhylogeneticModel)
varnames(PM::GroupBasedPhylogeneticModel)Return the symbolic name for the probability coordinates if PM is a PhylogeneticModel and for the Fourier corrdinates if PM is a GroupBasedPhylogeneticModel.
Examples
julia> tree = phylogenetic_tree(graph_from_edges(Directed,[[4,1], [4,2], [4,3]]));
julia> PM = general_markov_model(tree);
julia> varnames(PM)
"p"
julia> varnames(jukes_cantor_model(tree))
"q"This function is part of the experimental code in Oscar. Please read here for more details.
Phylogenetic rings
Two main rings are associated with any phylogenetic model: the parameter ring, which contains the unobserved parameters of the evolutionary process (e.g., transition probabilities), and the model ring, which represents the observed data (e.g., joint probability distribution at the leaves).
Parameter ring
The parameter ring's generators correspond to the symbolic parameters of the model. For a PhylogeneticModel, these are the entries of the transition matrices and root distribution. For a GroupBasedPhylogeneticModel, they are the Fourier parameters.
parameter_ring — Method
parameter_ring(PM::PhylogeneticModel; cached=false, sorted_edges::Union{Vector{Edge}, Nothing}=nothing)Create the polynomial ring for the parameters of a phylogenetic model PhylogeneticModel{GT, L, M, T} where GT <: PhylogeneticTree or GT <: PhylogeneticNetwork.
The optional sorted_edges argument is a vector of Edges that specifies a precise ordering for the edges when creating the polynomial variables for the model parameters. If sorted_edges is not provided (nothing), a default order is applied: the pendant edges corresponding to leaves 1 through n are listed first, followed by the internal edges, which are sorted based on their descendant nodes.
If M <: VarName, returns a tuple containing:
- The polynomial ring.
- A dictionary mapping parameter variables and edges to the corresponding ring generators.
- The root distribution vector.
If M <: MPolyRingElem, then returns a tuple containing:
- The polynomial ring.
- A dictionary mapping each edge to a homomorphism (
Oscar.MPolyAnyMap) from the original transition matrix ring to the new parameter ring. - The root distribution vector.
If GT <: PhylogeneticNetwork it additionally returns:
- A dictionary mapping hybrid edges to the corresponding hybrid parameter in the ring.
This function is part of the experimental code in Oscar. Please read here for more details.
parameter_ring — Method
parameter_ring(PM::GroupBasedPhylogeneticModel; cached=false, sorted_edges::Union{Vector{Edge}, Nothing}=nothing)Create the polynomial ring for the Fourier parameters of the model.
The optional sorted_edges argument is a vector of Edges that specifies a precise ordering for the edges when creating the polynomial variables for the model parameters. If sorted_edges is not provided (nothing), a default order is applied: the pendant edges corresponding to leaves 1 through n are listed first, followed by the internal edges, which are sorted based on their descendant nodes.
Returns a tuple containing:
- The polynomial ring.
- A dictionary mapping parameter variables and edges to the corresponding ring generators.
If GT <: PhylogeneticNetwork it additionally returns:
- A dictionary mapping hybrid edges to the corresponding hybrid parameter in the ring.
This function is part of the experimental code in Oscar. Please read here for more details.
To access specific parameters (entries of the treansition matrices, fourier parameters, entry of the root distriution or paramters associated to hybrid nodes in phylo networks) you can use the following functions:
entry_root_distribution — Function
entry_root_distribution(PM::PhylogeneticModel, i::Int)
entry_root_distribution(PM::GroupBasedPhylogeneticModel, i::Int)Return the i-th entry of the root distribution parameter for the phylogenetic model associated to PM.
This function is part of the experimental code in Oscar. Please read here for more details.
entry_transition_matrix — Function
entry_transition_matrix(PM::PhylogeneticModel, i::Int, j::Int, e::Edge)
entry_transition_matrix(PM::PhylogeneticModel, i::Int, j::Int, u::Int, v::Int)Return the (i, j)-th entry of the transition matrix associated with edge e (or equivalently Edge(u,v)) in the phylogenetic model PM.
This function is part of the experimental code in Oscar. Please read here for more details.
entry_transition_matrix(PM::GroupBasedPhylogeneticModel, i::Int, j::Int, e::Edge)
entry_transition_matrix(PM::GroupBasedPhylogeneticModel, i::Int, j::Int, u::Int, v::Int)Return the (i, j)-th entry of the transition matrix associated with edge e (or equivalently Edge(u,v)) of the underlying PhylogeneticModel contained within a GroupBasedPhylogeneticModel.
This function is part of the experimental code in Oscar. Please read here for more details.
entry_fourier_parameter — Function
entry_fourier_parameter(PM::GroupBasedPhylogeneticModel, i::Int, e::Edge)Return the i-th Fourier parameter associated with edge e (or equivalently Edge(u,v)) in a group-based phylogenetic model.
This function is part of the experimental code in Oscar. Please read here for more details.
entry_hybrid_parameter — Function
entry_hybrid_parameter(PM::Union{GroupBasedPhylogeneticModel{<: PhylogeneticNetwork}, PhylogeneticModel{<: PhylogeneticNetwork}}, e::Edge)
entry_hybrid_parameter(PM::Union{GroupBasedPhylogeneticModel{<: PhylogeneticNetwork}, PhylogeneticModel{<: PhylogeneticNetwork}}, u::Int, v::Int)Return the parameter associated with a hybrid edge e (or edge Edge(u,v)) in a phylogenetic model or groupd-based model PM.
This function is part of the experimental code in Oscar. Please read here for more details.
For example
julia> tree = graph_from_edges(Directed,[[4,1], [4,2], [4,3]]);
julia> PM = kimura2_model(tree)
Group-based phylogenetic model on a tree with 3 leaves and 3 edges
with root distribution [1//4, 1//4, 1//4, 1//4],
transition matrices of the form
[:a :b :c :b;
:b :a :b :c;
:c :b :a :b;
:b :c :b :a]
and fourier parameters of the form [:x, :y, :z, :z].
julia> S, _ = parameter_ring(PM);
julia> S
Multivariate polynomial ring in 9 variables x[1], x[2], x[3], y[1], ..., z[3]
over rational field
julia> entry_root_distribution(PM, 1)
1//4
julia> entry_transition_matrix(PM, 3, 3, Edge(4, 1))
a[1]
julia> entry_fourier_parameter(PM, 4, Edge(4, 1))
z[1]Model Ring and Equivalence Classes
The model ring's variables represent the joint probability distribution of states at the leaves of the tree. A key concept is that different leaf configurations can yield the same probability polynomial under the model's parametrization. These sets of configurations are called equivalence classes.
The fullmodelring contains a generator for every possible configuration of states at the leaves. A more compact and often more useful ring is the model_ring, which has one generator for each equivalence class.
full_model_ring — Method
full_model_ring(PM::Union{PhylogeneticModel, GroupBasedPhylogeneticModel}; cached=false)Creates the full model ring in probability coordinates for a PhylogeneticModel. Equivalently, creates the full model ring in the Fourier coordinates for a GroupBasedPhylogeneticModel.
The output consist on a tuple containing a IndexedRing and a Dict mapping a tuple of states at the leaves to the corresponding coordinate in the ring. This ring has a generator for every possible state configuration at the leaves.
Examples
julia> tree = phylogenetic_tree(graph_from_edges(Directed,[[4,1], [4,2], [4,3]]));
julia> R, p = full_model_ring(general_markov_model(tree));
julia> p[1,2,4]
p[1,2,4]This function is part of the experimental code in Oscar. Please read here for more details.
equivalent_classes — Method
equivalent_classes(PM::Union{PhylogeneticModel, GroupBasedPhylogeneticModel})Calculates the equivalence classes of the model's parametrizations. Equivalent polynomials are those that are identical under the chosen parametrization. This is a crucial step for dimension reduction in the model's algebraic variety.
Examples
julia> tree = phylogenetic_tree(graph_from_edges(Directed,[[4,1], [4,2], [4,3]]));
julia> PM = jukes_cantor_model(tree);
julia> R, q = full_model_ring(PM);
julia> class = equivalent_classes(PM);
julia> class[1,2,2]
3-element Vector{MPolyRingElem}:
q[1,2,2]
q[1,3,3]
q[1,4,4]This function is part of the experimental code in Oscar. Please read here for more details.
model_ring — Method
model_ring(PM::Union{PhylogeneticModel, GroupBasedPhylogeneticModel}; cached=false)Creates a reduced model ring based on the equivalence classes of the parametrizations. This ring has a generator for each unique polynomial (i.e., each equivalence class).
Examples
julia> tree = phylogenetic_tree(graph_from_edges(Directed,[[4,1], [4,2], [4,3]]));
julia> R, q = model_ring(cavender_farris_neyman_model(tree));
julia> q[1,1,1]
q[1,1,1]This function is part of the experimental code in Oscar. Please read here for more details.
Parametrizations
A parametrization is a ring homomorphism from the model ring (representing observed probabilities) to the parameter ring (representing model parameters). This map explicitly gives the polynomial formula for each observable probability in terms of the model parameters.
There are two versions: full_parametrization maps from the full_model_ring, while parametrization maps from the reduced model_ring.
full_parametrization — Function
full_parametrization(PM::PhylogeneticModel{<:PhylogeneticTree})
full_parametrization(PM::PhylogeneticModel{<:PhylogeneticNetwork})Constructs the parametrization map from the full model ring in probability coordinates (with generators for all leaf probability configurations) to the parameter ring (with generators from the transition matrices and the root distribution).
This function is part of the experimental code in Oscar. Please read here for more details.
full_parametrization(PM::GroupBasedPhylogeneticModel{<:PhylogeneticTree})
full_parametrization(PM::GroupBasedPhylogeneticModel{<:PhylogeneticNetwork})Constructs the parametrization map from the full model ring in Fourier coordinates (with generators for all leaf probability configurations) to the parameter ring (with generators from the Fourier parameters).
This function is part of the experimental code in Oscar. Please read here for more details.
parametrization — Method
parametrization(PM::Union{PhylogeneticModel, GroupBasedPhylogeneticModel})Constructs the parametrization map from the model ring to the parameter ring.
The rings used depend on the type of PM:
PhylogeneticModel: Map from probability coordinates to parameters of the transition matrices and root distribution.GroupBasedPhylogeneticModel: Map from Fourier coordinates to Fourier parameters.
This function is part of the experimental code in Oscar. Please read here for more details.
Affine Parametrization: Sometimes one wants to work on an affine variety by imposing constraints, such as the rows of the transition matrices summing to one. affine_parametrization provides these constrained maps.
affine_parametrization — Function
affine_parametrization(PM::PhylogeneticModel)
full_affine_parametrization(PM::PhylogeneticModel)Constructs an affine parametrization (reduced or full) of the model by imposing the transition matrices rows sum to one and the root distribution sums to one.
Examples
julia> tree = phylogenetic_tree(graph_from_edges(Directed,[[4,1], [4,2], [4,3]]));
julia> PM = phylogenetic_model(jukes_cantor_model(tree));
julia> affine_parametrization(PM)
Ring homomorphism
from multivariate polynomial ring in 5 variables over QQ
to multivariate polynomial ring in 6 variables over QQ
defined by
p[1,2,3] -> -2*b[1]*b[2]*b[3] + 1//4*b[1]*b[2] + 1//4*b[1]*b[3] + 1//4*b[2]*b[3]
p[1,2,2] -> 2*b[1]*b[2]*b[3] - 3//4*b[1]*b[2] - 3//4*b[1]*b[3] + 1//4*b[1] + 1//4*b[2]*b[3]
p[1,2,1] -> 2*b[1]*b[2]*b[3] - 3//4*b[1]*b[2] + 1//4*b[1]*b[3] - 3//4*b[2]*b[3] + 1//4*b[2]
p[1,1,2] -> 2*b[1]*b[2]*b[3] + 1//4*b[1]*b[2] - 3//4*b[1]*b[3] - 3//4*b[2]*b[3] + 1//4*b[3]
p[1,1,1] -> -6*b[1]*b[2]*b[3] + 9//4*b[1]*b[2] + 9//4*b[1]*b[3] - 3//4*b[1] + 9//4*b[2]*b[3] - 3//4*b[2] - 3//4*b[3] + 1//4This function is part of the experimental code in Oscar. Please read here for more details.
affine_parametrization(PM::GroupBasedPhylogeneticModel)Constructs an affine parametrization of the group-based model by setting the first Fourier parameters to 1.
Examples
julia> tree = phylogenetic_tree(graph_from_edges(Directed,[[4,1], [4,2], [4,3]]));
julia> PM = jukes_cantor_model(tree);
julia> affine_parametrization(PM)
Ring homomorphism
from multivariate polynomial ring in 5 variables over QQ
to multivariate polynomial ring in 6 variables over QQ
defined by
q[2,3,4] -> y[1]*y[2]*y[3]
q[2,2,1] -> y[1]*y[2]
q[2,1,2] -> y[1]*y[3]
q[1,2,2] -> y[2]*y[3]
q[1,1,1] -> 1This function is part of the experimental code in Oscar. Please read here for more details.
Example of a full workflow:
julia> tree = graph_from_edges(Directed,[[4,1], [4,2], [4,3]]);
julia> JC = jukes_cantor_model(tree);
julia> pmJC = phylogenetic_model(JC);
julia> Rp, p = model_ring(pmJC);
julia> Rq, q = model_ring(JC);
julia> f = parametrization(pmJC);
julia> f(p[(1, 1, 1)]) # Probability of observing state 1 at all leaves
1//4*a[1]*a[2]*a[3] + 3//4*b[1]*b[2]*b[3]
julia> g = parametrization(JC);
julia> g(q[(1, 1, 1)])
x[1]*x[2]*x[3]
julia> vanishing_ideal(pmJC)
Ideal generated by
2*p[1,2,3]^3 - p[1,2,3]^2*p[1,2,2] - p[1,2,3]^2*p[1,2,1] - p[1,2,3]^2*p[1,1,2] + p[1,2,3]^2*p[1,1,1] - p[1,2,3]*p[1,2,2]^2 - p[1,2,3]*p[1,2,1]^2 - p[1,2,3]*p[1,1,2]^2 + p[1,2,3]*p[1,1,1]^2 + p[1,2,2]^2*p[1,2,1] + p[1,2,2]^2*p[1,1,2] + p[1,2,2]*p[1,2,1]^2 - p[1,2,2]*p[1,2,1]*p[1,1,2] - p[1,2,2]*p[1,2,1]*p[1,1,1] + p[1,2,2]*p[1,1,2]^2 - p[1,2,2]*p[1,1,2]*p[1,1,1] + p[1,2,1]^2*p[1,1,2] + p[1,2,1]*p[1,1,2]^2 - p[1,2,1]*p[1,1,2]*p[1,1,1]
julia> vanishing_ideal(JC)
Ideal generated by
-q[2,3,4]^2*q[1,1,1] + q[2,2,1]*q[2,1,2]*q[1,2,2]The vanishing ideal can be computed using different algorithms
julia> tree = graph_from_edges(Directed,[[4,1], [4,2], [4,3]]);
julia> JC = jukes_cantor_model(tree);
julia> vanishing_ideal(JC; algorithm = :kernel)
Ideal generated by
-q[2,3,4]^2*q[1,1,1] + q[2,2,1]*q[2,1,2]*q[1,2,2]
julia> vanishing_ideal(JC; algorithm = :eliminate)
Ideal generated by
-q[2,3,4]^2*q[1,1,1] + q[2,2,1]*q[2,1,2]*q[1,2,2]
julia> vanishing_ideal(JC; algorithm = :f4)
Ideal generated by
-q[2,3,4]^2*q[1,1,1] + q[2,2,1]*q[2,1,2]*q[1,2,2]
julia> vanishing_ideal(JC; algorithm = :markov)
Ideal generated by
q[2,3,4]^2*q[1,1,1] - q[2,2,1]*q[2,1,2]*q[1,2,2]Coordinate Change for Group-Based Models
For GroupBasedPhylogeneticModel objects, this module provides functions to convert between the standard probability coordinates and the Fourier coordinates. This change of basis often simplifies the phylogenetic invariants and the geometry of the model.
The functions fourier_transform and inverse_fourier_transform compute the transformation matrices, while coordinate_change and inverse_coordinate_change return the actual ring homomorphisms.
fourier_transform — Method
fourier_transform(PM::GroupBasedPhylogeneticModel)Computes specialized Fourier transform from the matrix that represents the full Fourier transform. This matrix transforms the reduced probability coordinates (corresponding to the equivalent classes) to the reduced Fourier coordinates.
This function is part of the experimental code in Oscar. Please read here for more details.
inverse_fourier_transform — Method
inverse_fourier_transform(PM::GroupBasedPhylogeneticModel)Computes the matrix that transforms Fourier coordinates back to probability coordinates.
This function is part of the experimental code in Oscar. Please read here for more details.
coordinate_change — Method
coordinate_change(PM::GroupBasedPhylogeneticModel)Constructs the map from probability coordinates to Fourier coordinates. This is a homomorphism between the model rings.
This function is part of the experimental code in Oscar. Please read here for more details.
inverse_coordinate_change — Method
inverse_coordinate_change(PM::GroupBasedPhylogeneticModel)Constructs the map from Fourier coordinates to probability coordinates. This is the inverse of coordinate_change.
This function is part of the experimental code in Oscar. Please read here for more details.
julia> tree = graph_from_edges(Directed,[[4,1], [4,2], [4,3]]);
julia> CFN = cavender_farris_neyman_model(tree);
julia> cc = coordinate_change(CFN);
julia> cc_inv = inverse_coordinate_change(CFN);
julia> Rq, q = model_ring(CFN);
julia> Rp, p = model_ring(phylogenetic_model(CFN));
julia> cc(q[(1, 1, 1)])
p[1,2,2] + p[1,2,1] + p[1,1,2] + p[1,1,1]
julia> cc_inv(cc(q[(1, 1, 1)]))
q[1,1,1]Phylogenetic Networks
The PhylogeneticNetwork structure provides a way to represent and manipulate phylogenetic networks, which are directed acyclic graphs (DAGs) that generalizes phylogenetic trees by allowing hybrid nodes — nodes with multiple incoming edges — that represent genetic mixing events between lineages.
PhylogeneticNetwork — Type
PhylogeneticNetwork{N, L} <: AbstractGraph{Directed}A data structure representing a phylogenetic network.
A phylogenetic network is a directed acyclic graph where leaves (nodes with out-degree 0) represent observed species, and internal nodes represent ancestral species. Hybrid nodes (nodes with in-degree greater than 1) represent reticulation events like hybridization.
Type Parameters:
N: The number of hybrid nodes in the network.L: The level of the phylogenetic network.
Fields:
graph::Graph{Directed}: The underlying directed graph representing the network topology.hybrids::Dict{Int, Vector{Edge}}: A dictionary mapping each hybrid vertex ID to its incoming hybrid edges.
This function is part of the experimental code in Oscar. Please read here for more details.
phylogenetic_network — Function
phylogenetic_network(G::Graph{Directed})A constructor for creating a PhylogeneticNetwork from a directed graph G. It automatically detects hybrid nodes and edges.
Example
julia> G1 = graph_from_edges(Directed,[[5,6], [5,4], [6,4], [5,2], [6,3], [4,1]]);
julia> phylogenetic_network(G1)
Level-1 phylogenetic network with hybrid nodes {4} and edges
(4, 1)(5, 2)(5, 4)(5, 6)(6, 3)(6, 4)
julia> G2 = graph_from_edges(Directed,[[7,2], [8,7], [7,6], [8,6], [9,8], [10,9], [11,10], [11,12], [12,10], [7,2], [6,1], [9,3], [11,4], [12,5]]);
julia> phylogenetic_network(G2)
Level-1 phylogenetic network with hybrid nodes {6, 10} and edges
(6, 1)(7, 2)(7, 6)(8, 6)(8, 7)(9, 3)(9, 8)(10, 9)(11, 4)(11, 10)(11, 12)(12, 5)(12, 10)
julia> G3 = graph_from_edges(Directed,[[6,2], [6,7], [6,5], [7,5], [5,4], [7,8], [8,4], [4,1], [8,3]]);
julia> phylogenetic_network(G3)
Level-2 phylogenetic network with hybrid nodes {4, 5} and edges
(4, 1)(5, 4)(6, 2)(6, 5)(6, 7)(7, 5)(7, 8)(8, 3)(8, 4)This function is part of the experimental code in Oscar. Please read here for more details.
You can access specific properties of a PhylogeneticNetwork, like its level and hybrid nodes, using these functions:
level — Method
level(::PhylogeneticNetwork{N, L}) where {N, L}Return the level L of the phylogenetic network.
Example
julia> N = Oscar.phylogenetic_network(graph_from_edges(Directed,[[5,6], [5,4], [6,4], [5,2], [6,3], [4,1]]));
julia> level(N)
1This function is part of the experimental code in Oscar. Please read here for more details.
n_hybrid — Function
n_hybrid(::PhylogeneticNetwork{N, L}) where {N, L}Return the number of hybrid nodes N in the phylogenetic network.
Example
julia> N = phylogenetic_network(graph_from_edges(Directed,[[5,6], [5,4], [6,4], [5,2], [6,3], [4,1]]));
julia> n_hybrid(N)
1 This function is part of the experimental code in Oscar. Please read here for more details.
hybrids — Function
hybrids(N::PhylogeneticNetwork)Return the dictionary mapping hybrid vertex IDs to their incoming hybrid edges.
Examples
julia> N = phylogenetic_network(graph_from_edges(Directed,[[5,6], [5,4], [6,4], [5,2], [6,3], [4,1]]));
julia> hybrids(N)
Dict{Int64, Vector{Edge}} with 1 entry:
4 => [Edge(5, 4), Edge(6, 4)]This function is part of the experimental code in Oscar. Please read here for more details.
hybrid_vertices — Function
hybrid_vertices(N::PhylogeneticNetwork)Return a list of the hybrid vertices in the phylogenetic network.
Examples
julia> N = phylogenetic_network(graph_from_edges(Directed,[[5,6], [5,4], [6,4], [5,2], [6,3], [4,1]]));
julia> hybrid_vertices(N)
1-element Vector{Int64}:
4This function is part of the experimental code in Oscar. Please read here for more details.
hybrid_edges — Function
hybrid_edges(N::PhylogeneticNetwork)Return a list of all hybrid edges in the phylogenetic network.
Examples
julia> N = phylogenetic_network(graph_from_edges(Directed,[[5,6], [5,4], [6,4], [5,2], [6,3], [4,1]]));
julia> hybrid_edges(N)
1-element Vector{Vector{Edge}}:
[Edge(5, 4), Edge(6, 4)]This function is part of the experimental code in Oscar. Please read here for more details.
General graph properties can be accessed with the following functions:
graph — Method
graph(N::PhylogeneticNetwork)Return the underlying Graph{Directed} of the phylogenetic network.
Examples
julia> N = phylogenetic_network(graph_from_edges(Directed,[[5,6], [5,4], [6,4], [5,2], [6,3], [4,1]]));
julia> graph(N)
Directed graph with 6 nodes and the following edges:
(4, 1)(5, 2)(5, 4)(5, 6)(6, 3)(6, 4)This function is part of the experimental code in Oscar. Please read here for more details.
n_edges — Method
n_edges(N::PhylogeneticNetwork)Return the total number of edges in the phylogenetic network.
Examples
julia> N = phylogenetic_network(graph_from_edges(Directed,[[5,6], [5,4], [6,4], [5,2], [6,3], [4,1]]));
julia> n_edges(N)
6This function is part of the experimental code in Oscar. Please read here for more details.
This function is part of the experimental code in Oscar. Please read here for more details.
This function is part of the experimental code in Oscar. Please read here for more details.
edges — Method
edges(N::PhylogeneticNetwork)Return an iterator over the edges of the phylogenetic network.
Examples
julia> N = phylogenetic_network(graph_from_edges(Directed,[[5,6], [5,4], [6,4], [5,2], [6,3], [4,1]]));
julia> collect(edges(N))
6-element Vector{Edge}:
Edge(4, 1)
Edge(5, 2)
Edge(5, 4)
Edge(5, 6)
Edge(6, 3)
Edge(6, 4)This function is part of the experimental code in Oscar. Please read here for more details.
leaves — Method
leaves(N::PhylogeneticNetwork)Return the leaf nodes of the phylogenetic network.
Examples
julia> N = phylogenetic_network(graph_from_edges(Directed,[[5,6], [5,4], [6,4], [5,2], [6,3], [4,1]]));
julia> leaves(N)
3-element Vector{Int64}:
1
2
3This function is part of the experimental code in Oscar. Please read here for more details.
All model constructors, including PhylogeneticModel, GroupBasedPhylogeneticModel, and the classic pre-defined models, can be applied to phylogenetic networks in the same way they are applied to trees. The resulting model object will be defined on the network, incorporating its specific topology. Consequently, all network property functions (e.g., graph, level, hybrids) can be called directly on the model object.
Example: CFN model on a Level-1 Network
The following example demonstrates this by creating a Cavender-Farris-Neyman model on a level-1 network. It shows how to inspect the model's underlying graph and access the algebraic parameters. We use the method components_of_kernel(d::Int, phi::MPolyAnyMap) with $d=3$ to compute all minimal generators of the kernel of the polynomial map $\phi$ with total degree at most $3$ using the main algorithm of [CH26].
julia> G = graph_from_edges(Directed,[[7,6], [7,8], [6,5], [8,5], [5,1], [6,2], [7,3], [8,4]]);
julia> PM = cavender_farris_neyman_model(G)
Group-based phylogenetic model on a level-1 network with 1 hybrid node, 4 leaves
and 8 edges with root distribution [1//2, 1//2],
transition matrices of the form
[:a :b;
:b :a]
and fourier parameters of the form [:x, :y].
julia> graph(PM)
Level-1 phylogenetic network with hybrid nodes {5} and edges
(5, 1)(6, 2)(6, 5)(7, 3)(7, 6)(7, 8)(8, 4)(8, 5)
julia> S, _ = parameter_ring(PM);
julia> S
Multivariate polynomial ring in 18 variables l[1, 1], l[1, 2], x[1], x[2], ..., y[8]
over rational field
julia> entry_fourier_parameter(PM, 1, Edge(6, 5))
x[5]
julia> entry_hybrid_parameter(PM, Edge(6, 5))
l[1, 1]
julia> R, q = full_model_ring(PM);
julia> f = full_parametrization(PM);
julia> f(q[1,1,1,1])
l[1, 1]*x[1]*x[2]*x[3]*x[4]*x[5]*x[6]*x[7] + l[1, 2]*x[1]*x[2]*x[3]*x[4]*x[6]*x[7]*x[8]
julia> phi = parametrization(PM)
Ring homomorphism
from multivariate polynomial ring in 8 variables over QQ
to multivariate polynomial ring in 18 variables over QQ
defined by
q[2,2,2,2] -> l[1, 1]*x[6]*y[1]*y[2]*y[3]*y[4]*y[5]*y[7] + l[1, 2]*x[7]*y[1]*y[2]*y[3]*y[4]*y[6]*y[8]
q[2,2,1,1] -> l[1, 1]*x[3]*x[4]*x[6]*x[7]*y[1]*y[2]*y[5] + l[1, 2]*x[3]*x[4]*y[1]*y[2]*y[6]*y[7]*y[8]
q[2,1,2,1] -> l[1, 1]*x[2]*x[4]*x[7]*y[1]*y[3]*y[5]*y[6] + l[1, 2]*x[2]*x[4]*x[6]*y[1]*y[3]*y[7]*y[8]
q[2,1,1,2] -> l[1, 1]*x[2]*x[3]*y[1]*y[4]*y[5]*y[6]*y[7] + l[1, 2]*x[2]*x[3]*x[6]*x[7]*y[1]*y[4]*y[8]
q[1,2,2,1] -> l[1, 1]*x[1]*x[4]*x[5]*x[7]*y[2]*y[3]*y[6] + l[1, 2]*x[1]*x[4]*x[7]*x[8]*y[2]*y[3]*y[6]
q[1,2,1,2] -> l[1, 1]*x[1]*x[3]*x[5]*y[2]*y[4]*y[6]*y[7] + l[1, 2]*x[1]*x[3]*x[8]*y[2]*y[4]*y[6]*y[7]
q[1,1,2,2] -> l[1, 1]*x[1]*x[2]*x[5]*x[6]*y[3]*y[4]*y[7] + l[1, 2]*x[1]*x[2]*x[6]*x[8]*y[3]*y[4]*y[7]
q[1,1,1,1] -> l[1, 1]*x[1]*x[2]*x[3]*x[4]*x[5]*x[6]*x[7] + l[1, 2]*x[1]*x[2]*x[3]*x[4]*x[6]*x[7]*x[8]
julia> components = components_of_kernel(3, phi);
julia> length(components)
1
julia> coordinate_change(PM)
Ring homomorphism
from multivariate polynomial ring in 8 variables over QQ
to multivariate polynomial ring in 8 variables over QQ
defined by
q[2,2,2,2] -> -p[1,2,2,2] + p[1,2,2,1] + p[1,2,1,2] - p[1,2,1,1] + p[1,1,2,2] - p[1,1,2,1] - p[1,1,1,2] + p[1,1,1,1]
q[2,2,1,1] -> -p[1,2,2,2] - p[1,2,2,1] - p[1,2,1,2] - p[1,2,1,1] + p[1,1,2,2] + p[1,1,2,1] + p[1,1,1,2] + p[1,1,1,1]
q[2,1,2,1] -> -p[1,2,2,2] - p[1,2,2,1] + p[1,2,1,2] + p[1,2,1,1] - p[1,1,2,2] - p[1,1,2,1] + p[1,1,1,2] + p[1,1,1,1]
q[2,1,1,2] -> -p[1,2,2,2] + p[1,2,2,1] - p[1,2,1,2] + p[1,2,1,1] - p[1,1,2,2] + p[1,1,2,1] - p[1,1,1,2] + p[1,1,1,1]
q[1,2,2,1] -> p[1,2,2,2] + p[1,2,2,1] - p[1,2,1,2] - p[1,2,1,1] - p[1,1,2,2] - p[1,1,2,1] + p[1,1,1,2] + p[1,1,1,1]
q[1,2,1,2] -> p[1,2,2,2] - p[1,2,2,1] + p[1,2,1,2] - p[1,2,1,1] - p[1,1,2,2] + p[1,1,2,1] - p[1,1,1,2] + p[1,1,1,1]
q[1,1,2,2] -> p[1,2,2,2] - p[1,2,2,1] - p[1,2,1,2] + p[1,2,1,1] + p[1,1,2,2] - p[1,1,2,1] - p[1,1,1,2] + p[1,1,1,1]
q[1,1,1,1] -> p[1,2,2,2] + p[1,2,2,1] + p[1,2,1,2] + p[1,2,1,1] + p[1,1,2,2] + p[1,1,2,1] + p[1,1,1,2] + p[1,1,1,1]