 Research
 Open access
 Published:
Unifying duplication episode clustering and genespecies mapping inference
Algorithms for Molecular Biology volume 19, Article number: 7 (2024)
Abstract
We present a novel problem, called MetaEC, which aims to infer genespecies assignments in a collection of partially leaflabeled gene trees labels by minimizing the size of duplication episode clustering (EC). This problem is particularly relevant in metagenomics, where incomplete data often poses a challenge in the accurate reconstruction of gene histories. To solve MetaEC, we propose a polynomial time dynamic programming (DP) formulation that verifies the existence of a set of duplication episodes from a predefined set of episode candidates. In addition, we design a method to infer distributions of genespecies mappings. We then demonstrate how to use DP to design an algorithm that solves MetaEC. Although the algorithm is exponential in the worst case, we introduce a heuristic modification of the algorithm that provides a solution with the knowledge that it is exact. To evaluate our method, we perform two computational experiments on simulated and empirical data containing whole genome duplication events, showing that our algorithm is able to accurately infer the corresponding events.
Introduction
In the field of computational biology, the use of gene families and the reconciliation model has become increasingly popular for studying the evolution of diverse organisms. These tools have facilitated the development of new algorithms and computational methods capable of handling large and complex datasets and exploring various types of evolutionary events. These events range from simple macroevolutionary processes such as gene duplications, gene losses, and horizontal gene transfers to more complex ones such as genomic duplications, speciations, and hybridizations. The reconciliation model has enabled researchers to reconstruct evolutionary histories by reconciling gene trees with species trees and identifying the evolutionary events that have led to the observed patterns of gene evolution. In the context of metagenomics, the reconciliation model has also been used to detect missing genespecies assignments using polynomial time algorithms. These developments have led to a better understanding of the evolutionary processes.
A classical reconciliation model [1, 2] defines a mapping from every node from a gene tree into a node in the species tree and determines if such a node is related to a speciation or can be classified as a single gene duplication [3]. In result, an embedding of the set of gene trees into a species tree can be interpreted as a joint evolutionary scenario [4]. The classical least common ancestor (LCA) mapping minimizes the number of single gene duplications and losses for one gene tree and the species tree [4].
The whole genome duplication (WGD) phenomenon incorporates additional copies of a complete genome into the original genetic material, thus creating an opportunity to introduce novel evolutionary traits [5, 6]. From a macro perspective, this phenomenon played a crucial role in the divergence and formation of species and shaped the evolution of almost all major lineages of life. In particular, many WGDs were uncovered in the evolutionary histories of plants, especially crops. WDGs potentially enabled the successful domestication of plants [7] and are important in the fight against famine [8]. Many traces and evidence of whole genome duplications can be found in the genomes of yeast and other fungal species [5, 9]. From the perspective of single cell evolution studies, WGDs are prevalent in cancer progression [10] and can lead us to the prognosis of advanced cancer stages [11] or the creation of strategies for targeted therapy [12].
Guigó et al. [13] proposed the first approach for detecting multiple gene duplication episodes from a collection of rooted gene trees. They designed a heuristic that aggregates single gene duplication events into a large gene duplication, given a collection of rooted gene trees and a rooted species tree. This approach was formalized and improved by Page and Cotton [14], who defined the problem of episode clustering (\(\textsf{EC}\)) as the task of identifying the minimal number of locations in the species tree where all duplications from the input gene trees can be placed. Fellows [15] applied this model in the context of the supertree problem. Polynomialtime solutions for two types of multiple gene duplication problems episode clustering and a more general variant of clustering called minimum episodes (ME) were proposed in [16, 17]. Luo et al. [18] proposed linear time and space algorithms, partially based on [19], for these problems. [20] introduced a unified approach by proposing a concept of interval models with a linear time and space solution to a broad class of clustering problems including EC and ME. Alternative approaches include generalization to unrooted gene trees; however, such approaches are often computationally complex [21, 22]. Other approaches include variants of clustering rules that depend on the maximal number of duplication episodes placed in one path [23, 24]. A comprehensive analysis of various models is available in [20]. Furthermore, [25] proposes an integer linear programming formulation that simplifies the process of testing these models. Relevant computational complexity results on the ME problem are presented in [26].
Metagenomic studies provide valuable information for analyzing entire communities of organisms and revealing a complete picture of their functional and adaptive capacities crucial for ecology [27] or human health [28]. Genetic material isolated in such studies can be used to detect whole genome duplication events.
One of the steps in metagenomic analysis is called binning. The aim of this procedure is to assign sequenced DNA fragments to the appropriate taxonomic groups. The assignment of certain genes to species can be ambiguous due to the limitations of annotation methods. A precise and comprehensive gene tree topology is essential for the accurate identification of potential duplication sites. The absence or misplacement of duplications in gene trees can, in turn, result in incorrect outcomes of methods aimed at determining whole genome duplication events.
To tackle the challenge of missing genespecies assignments in evolutionary studies, previous research has introduced methods based on the reconciliation score using gene duplication and loss events [29, 30]. In a related work, Mykowiecka et al. [31] extended this model by including horizontal gene transfer to better analyze bacterial evolution and proposed polynomial time algorithms for these models. These approaches utilize tree reconciliation according to the classical scheme, in which the gene tree includes symbols representing sequences with unclear species assignment in addition to the known gene labels. The objective is to assign the unknown gene labels to their corresponding species in a gene tree while minimizing the total reconciliation score, which is typically a weighted sum of evolutionary events such as gene duplication, gene loss, and horizontal gene transfer.
Here, we present a novel problem called MetaEC, which aims to infer genespecies assignments in a collection of gene trees with unknown labels by minimizing the size of episode clustering. This problem is particularly relevant in metagenomics, where incomplete data often poses a challenge in the accurate reconstruction of gene histories. To solve MetaEC, we propose a dynamic programming (DP) algorithm that verifies the existence of a set of duplication episodes from a predefined set of episode candidates. We then demonstrate how to use DP to design an algorithm that solves MetaEC. In addition, we design an algorithm to infer distributions of genespecies mappings from the set of all optimal solutions inferred by DP. Although the algorithm is exponential in the worst case, we introduce a heuristic modification of the algorithm that provides a solution with the knowledge that it is exact. To evaluate our method, we perform two computational experiments on simulated and empirical data containing WGD events, showing that our algorithm can accurately infer the corresponding events. It is important to note that a direct comparison of our novel algorithm with alternative methods is not possible, as there is currently no existing approach for the simultaneous inference of duplications and genespecies mappings.
Basic definitions
Gene trees, species trees, and the duplication cost
We begin by recalling some basic definitions from graph theory. All trees in this article are rooted and binary, therefore we refer to them as trees. For a tree \(T=(V(T),E(T))\), by \(\textsf{root} (T)\) we denote the root, and by L(T) we denote the set of all leaves. Every nonleaf node will be called internal. A species tree is a tree whose leaves are called species. A gene tree over a species tree S is a tree with leaves labeled by the species from S. The set of all species present in a species tree or a gene tree T is denoted by \(\mathcal {L} (T)\). Note that for a species tree S, \(L(S)=\mathcal {L} (S)\). Also, for a gene tree G over S, \(\mathcal {L} (G) \subseteq L(S)\).
For nodes a and b, \(a \preceq b\) means that a and b are on the same path from the root, with b being closer to the root than a. We write \(a \prec b\) if \(a \preceq b\) and \(a \ne b\).
For a gene tree G over a species tree S, the least common ancestor (lca) mapping between G and S is a function \(\textsf{M} _{G} :V(G) \rightarrow V(S)\) defined as follows. If v is a leaf in G then \(\textsf{M} _G(v)\) is the label of v. When v is an internal node in T having two children a and b, then \(\textsf{M} _G(v)\) is the least common ancestor of \(\textsf{M} _G(a)\) and \(\textsf{M} _G(b)\) in S. An internal node \(g \in V(G)\) is called a duplication if \(\textsf{M} _G(g) = \textsf{M} _G(a)\) for a child a of g. The duplication cost, denoted by \(\textsf{D} (G,S)\), is the total number of duplications in G. Each nonduplication internal node of G we call a speciation. In the latter part of the article, a duplication g is called an sduplication if \(\textsf{M} _G(g)=s\). Similarly, we use the notation for an sspeciation and an sleaf.
An example of tree reconciliation and the lcamapping is depicted in the leftmost part of Fig. 1.
Episode clustering problems
Below we present a model of duplication episodes proposed in [20]. In short, this model admits all evolutionary scenarios using duplication and loss events with a minimal number of gene duplications.
Formally, the model of gene duplication episodes allows for relocating a gene duplication from its lcamapping node to one of its ancestors under some additional constraints required to preserve the biological soundness of the scenario. For a gene tree G over S, a mapping \(\textsf{F} _G :V(G) \rightarrow V(S)\) is called valid if the following conditions are satisfied:

\(\textsf{F} _G(a) \preceq \textsf{F} _G(b)\) if \(a \preceq b\) (time consistency),

\(\textsf{F} _G(a)=\textsf{M} _G(a)\) for any speciation node a (fixed speciations),

\(\textsf{F} _G(a) \succeq \textsf{M} _G(a)\) for any duplication node a (duplication can be raised),

\(\textsf{F} _G(a) \prec \textsf{M} _G(b)\) for any speciation node b such that \(a \prec b\) (fixed number of gene duplications).
Note that the model of valid mappings described above is more comprehensive than the model presented in [16].
Figure 1 provides an example of valid mappings that uniquely define an evolutionary scenario that can be represented as a tree with an additional decoration of nodes. For more information on the formal modeling of evolutionary scenarios, refer to [4].
We denote by \(\textsf{Dup} _T \subset V(T)\), the set of all duplication nodes in T. Let \(G_1,G_2,\dots ,G_n\) be a collection of rooted gene trees over a species tree S. Assume that, for every \(i \in \{1,2,\dots ,n\}\), \(\textsf{F} _i\) is a valid mapping between \(G_i\) and S. Every element \(s \in \bigcup _i \textsf{F} _i[\textsf{Dup} _{G_i}]\)^{Footnote 1} denotes the location of multiple gene duplication events in S. We will refer to such locations as duplication episodes or simply episodes. Later on, we may also use the term episode to refer to the set of duplications that are mapped into it.
Problem 1
(Episode Clustering, EC) Given a collection of rooted gene trees \(G_1,G_2,\dots ,G_n\) over a species tree S. Compute the minimal number of duplication episodes, denoted by \(\textsf{EC} (G_1,G_2,\dots , G_n, S)\), in the set of all valid mappings \(\textsf{F} _1,\textsf{F} _2,\dots ,\textsf{F} _n\) such that \(\textsf{F} _i :V(G_i) \rightarrow V(S)\).
This problem can be solved in linear time and space [18].
Genespecies mappings
We present the main problem for joint reconstruction of genespecies mappings and minimizing the set of episodes.
A partial gene tree is a rooted binary tree where each leaf is labeled by a species or has no label. Let G be a partial gene tree over S. By \(\Lambda _G :L(G) \rightarrow L(S)\) we denote the partial leaf labeling function such that \(\Lambda _G(g)\) is the label (species) of the leaf g in G if defined. Note that any gene tree is a partial gene tree with the leaf labeling being a total function. If a leaf in G has no label, we say that the label is unknown and write \(\Lambda _G(g)=\bot\). We say that a gene tree \(G^*\) over S extends a partial gene tree G over S if G and \(G^*\) are isomorphic as graphs (i.e., \(V(G)=V(G^*)\) and \(E(G)=E(G^*)\)), and, \(\Lambda _{G^*}\) is a total function that extends \(\Lambda _{G}\).
Inferring labelings by minimizing episodes
Now, we present the problem of the simultaneous reconstruction of leaf labelings and duplication episodes from collections of partial gene trees.
Problem 2
(MetaEC) Given a collection of partial gene trees \(G_1, G_2, \dots , G_k\) over a species tree S. Compute the minimum \(\textsf{EC} (G^*_1, G^*_2, \dots , G^*_k,S)\), denoted \(\textsf{EC} (G_1, G_2, \dots , G_k,S)\), in the set of all collections of gene trees \(G^*_1, G^*_2, \dots , G^*_k\) such that \(G^*_i\) extends \(G_i\), for each i.
For example, if \((a,(\bot ,(\bot ,\bot )))\) is a single gene tree with three unknown labels and (a, (b, c)) is a species tree, then the problem is to replace all occurrences of \(\bot\) by a, b or c such that the total number of duplication episodes is minimized. In this case, the optimal cost is 1, since at least one duplication is needed when the gene tree has four leaves and there are only three species.
Methods
We begin by solving a simpler problem in which we assume that the set of duplication episodes is constrained to a given set of species tree nodes. Next, we show how to solve MetaEC for a single gene tree. Section MetaEC in the general case presents the general solution.
Episode feasibility problem
We start with a related constrained problem. Given a partial gene tree, we are interested in the question of whether there is an extension of the partial gene tree such that the set of corresponding duplication episodes is contained in a given fixed set of episode candidates.
Problem 3
(Episode Feasibility) Given a partial gene tree G over a species tree S and \(X \subseteq V(S)\). Does there exist a gene tree \(G^*\) and a valid mapping \(\textsf{F} _{G^*}\) such that \(G^*\) extends G and \(\textsf{F} _{G^*}(\textsf{Dup} _{G^*}) \subseteq X\)?
If a partial gene tree G satisfies the above property, we call G Xfeasible with respect to a species tree S. If the context is clear, we omit the reference to S.
The solution to Episode Feasibility is a dynamic programming (DP) formulation expressed using Łukasiewicz’s ThreeValued Logic Ł\(_{3}\) [32] with three constants \(\textsf{True}\), \(\textsf{False}\), and \(\textsf{Unknown}\) (representing uncertainty) and ordered linearly: \(\textsf{False}< \textsf{Unknown} < \textsf{True}\). The logic has binary operators \(\vee\) (disjunction, \(\max\)), \(\wedge\) (conjunction, \(\min\)), and two unary operators \(\textsf{L}\) (certainty) and \(\textsf{M}\) (possibility). See the interpretation in Table 1.
For a node g of a gene tree G, by Gg we denote the subtree of G rooted at g. For any nonleaf node t in a tree, by \(t'\) and \(t''\), we denote the children of t. To simplify the notation, we assume that the set \(X \subseteq V(S)\) is fixed. Then, we have the following dynamic programming formulas that solve Episode Feasibility. Let \(g \in V(G)\) and \(s \in V(S)\).
where
In the next Lemma, we express properties satisfied by the above formulas.
For a partial gene tree G over S, and nodes \(g \in V(G)\) and \(s \in V(S)\), we say that a valid mapping \(\textsf{F} _T :V(Gg) \rightarrow V(Ss)\) is feasible for (g, s, X) if and only if T extends Gg and \(\textsf{F} _{T}(\textsf{Dup} _{T}) \subseteq X\). Feasible mappings represent episode scenarios that correspond to partial solutions to the instance of Episode Feasibility that forces duplications from Gg to be present in episodes from \(X \cap V(Ss)\).
We say a duplication g in a gene tree T is upper if the path from g to the root of T contains only duplications. The set of all upper duplications in a tree T is denoted by \(\textsf{UDup} _T\). We say that a valid mapping \(\textsf{F} _T :V(Gg) \rightarrow V(Ss)\) is weakly feasible for (g, s, X) if and only if there is no feasible mapping for (g, s, X), but there is T that extends Gg, T has at least one upper duplication d such that \(F_T(d) \notin X\) and \(\textsf{F} _{T}(\textsf{Dup} _T{\setminus }\textsf{UDup} _T) \subseteq X\). In contrast to feasible mappings, in weakly feasible mappings we constrain only nonupper duplications present in Gg. Here, the upper duplications are elements of episodes \(s \notin X\). This situation is modeled by \(\textsf{Unknown}\) value returned from \(\delta ^\downarrow (g,s)\) and \(\delta (g,s)\) calls, meaning that there is at least one duplication that needs to be assigned later (if possible) to an episode from \(X \setminus V(Ss)\), which eventually occurs at levels of recursion shallower than the level of (g, s).
Informally, the meaning of DP formulas can be understood as follows. Below, let T be an extension of the subtree of G rooted at a node g. The value of \(\delta (g, s)\) is \(\textsf{True}\) if there exists T where g is an sduplication and all duplications are assigned to the episodes from X (where \(s \in X\) as well). Similarly, the value of \(\sigma (g, s)\) is \(\textsf{True}\) if there exists T, where g represents an sspeciation or an sleaf node. Next, \(\delta (g, s)\) is \(\textsf{Unknown}\) if the condition for \(\delta (g, s) = \textsf{True}\) is not met. However, there still exists T, where g represents an sduplication, and all nonupper duplications from T are assigned to the episodes from X, while the upper duplications are assigned to episodes outside of X. It is important to note that in this case, s is not an element of X. Note that \(\sigma (g, s)\) cannot be \(\textsf{Unknown}\) since speciation nodes are fixed. Moving on, \(\delta ^\downarrow (g, s)\) is \(\textsf{True}\) if there exists T where all duplications are assigned to the episodes from X. Lastly, \(\delta ^\downarrow (g, s)\) is \(\textsf{Unknown}\) if the condition for \(\delta ^\downarrow (g, s) = \textsf{True}\) is not met. However, there still exists T, where all nonupper duplications from T are assigned to the episodes from X, while the upper duplications are assigned to episodes outside of X.
While \(\epsilon\) and \(\delta ^*\) should be treated as “local” in the main formulas (i.e., they should not form separate arrays in implementation), their properties can be formulated as follows. Generally, if \(\epsilon (g,s)\) is \(\textsf{True}\), then there is a tree T where g is mapped into s and all duplications are assigned to the episodes from X. Since \(\sigma (g,s)\) cannot be \(\textsf{Unknown}\), \(\epsilon (g, s)\) is \(\textsf{Unknown}\) only if \(\sigma (g,s)\) is \(\textsf{False}\) and \(\delta (g, s) = \textsf{Unknown}\). Again, here g is mapped to s. Now, \(\delta ^*(g,s)\) is \(\textsf{True}\) only if there is T where g is an sduplication, and all duplication nodes below g are assigned to episodes from X. Importantly, at least one of the children of g must be mapped to s, which is captured by \(\epsilon\). Furthermore, \(\delta ^*(g,s)\) resembles \(\delta (g,s)\), but the condition that duplications must be assigned to episodes from X only applies to the duplications (or upperduplications) below g if \(\delta ^*(g,s)\) is \(\textsf{True}\) (or \(\textsf{Unknown}\), respectively).
The following Lemma formalizes the conditions described above.
Lemma 1
Given a partial gene tree G over a species tree S and \(X \subseteq V(S)\). Let \(g \in V(G)\), \(s \in V(S)\). Then,

P1
\(\delta (g,s)=\textsf{True}\) if and only if there is a gene tree T and a feasible mapping \(\textsf{F} _T\) for (g, s, X) such that g is an sduplication in T.

P2
\(\delta (g,s)=\textsf{Unknown}\) if and only if there is a gene tree T and a weakly feasible mapping \(\textsf{F} _T\) for (g, s, X) such that g is an sduplication in T.

P3
\(\sigma (g,s)=\textsf{True}\) if and only if there a gene tree T and a feasible mapping \(\textsf{F} _T\) for (g, s, X) such that and g is an sspeciation or an sleaf in T.

P4
For any g and s, \(\sigma (g,s) \ne \textsf{Unknown}\).

P5
\(\delta ^\downarrow (g,s)\) is \(\textsf{True}\) if and only if there is a feasible mapping for (g, s, X).

P6
\(\delta ^\downarrow (g,s)\) is \(\textsf{Unknown}\) if and only if there is a weakly feasible mapping for (g, s, X).
Proof
The proof is by induction on the structure of G and S. The base of induction is when \(g \in L(G)\) and \(s \in L(S)\), for which all properties are easy to verify.
Inductive assumption: For every x, y such that \(g \succ x\), and \(s \succeq y\) or \(g \succeq x\) and \(s \succ y\), P1P6 are satisfied. Inductive hypothesis: For g and s, where at least one of g and s is not a leaf, P1P6 are satisfied.
Additional notation: since functions are relations, we identify a function \(f :A \rightarrow B\) with the sets of pairs \(\{ (x,f(x) :x \in A \}\). If \(\textsf{F} _{T'}:V(Gg') \rightarrow V(Ss)\) and \(\textsf{F} _{T''}:V(G g'')\rightarrow V(Ss)\) are valid mappings, then by \(\textsf{F} _{T'}\oplus \textsf{F} _{T''}\) we denote the mapping \(\textsf{F} _{T'}_{V(T') {\setminus } \textsf{UDup} _{T'}} \cup \textsf{F} _{T'}_{V(T'') {\setminus } \textsf{UDup} _{T''}} \cup \{(g,s)\} \cup (\textsf{UDup} _{T'} \cup \textsf{UDup} _{T''}) \times \{s\}\). Note that the resulting mapping is valid for a gene tree \((T',T'')\) that extends Gg.
We start with several properties.
(A1) if \(\delta ^\downarrow (g',s)=\textsf{Unknown}\), then \(s \notin X\). Assume that \(s \in X\). Then, by P6, there is a weakly feasible mapping \(\textsf{F} _{T'}\) for \((g',s,X)\). Now, let \(\textsf{F} _{T'}':= \textsf{F} _{T'}_{V(T') {\setminus } \textsf{UDup} _{T'}} \cup \textsf{UDup} _{T'} \times \{s\}\). It is not difficult to see that the mapping is feasible for \((g',s,X)\). A contradiction.
(A2) If \(\epsilon (g',s)=\textsf{Unknown}\), then there is a weakly feasible \(\textsf{F} _{T'}\) for \((g',s,X)\) and \(g'\) is an sduplication. Here, \(\sigma (g',s)\) cannot be \(\textsf{True}\), thus \(\sigma (g',s)=\textsf{False}\), by P3 and P4. Therefore, \(\delta (g',s)=\textsf{Unknown}\). The rest follows from P2.
We first prove properties P1P4 for \(\delta\) and \(\sigma\). Then, we show that P5 and P6 hold for g and s.
(P1, \(\Rightarrow\)): If \(\delta (g,s)=\textsf{True}\), then, from (1) g is internal, \(s \in X\) and \(\delta ^*(g,s) = \textsf{True}\). Then, w.l.o.g., for a child \(g'\) of g \(\epsilon (g',s) \wedge \delta ^\downarrow (g',s)=\textsf{True}\). Since, \(\delta (g',s) \vee \sigma (g',s)\) is \(\textsf{True}\), it follows from the inductive assumption for P1 and P3, that there is a feasible mapping \(\textsf{F} _{T'}\) for \((g',s,X)\) such that \(\textsf{M} _{T'}(g')=s\). For, the other child, we have \(\delta ^\downarrow (g'',s)=\textsf{True}\). From P5, there is a feasible mapping \(\textsf{F} _{T''}\) for \((g'',s,X)\). Now, let \(T=(T',T'')\). Then, the mapping \(\textsf{F} _{T'}\oplus \textsf{F} _{T''}\) is feasible for (g, s, X) and g is an sduplication in T.
(P1, \(\Leftarrow\)): Assume there is a feasible \(\textsf{F} _T\) for (g, s, X) such that \(\textsf{M} _T(g)=s\) and g is a duplication in T. Since \(g \in \textsf{UDup} _T\), we have that \(s \in X\) and only (1) is satisfied. Thus, \(\delta (g,s)=\delta ^*(g,s)\). W.l.o.g. we may assume that \(M_{T}(g')=s\) (recall that g is a duplication). Then, \(\textsf{F} _{T}_{V(Tg')}\) is feasible for \((g',s,X)\). By the inductive assumption for P1 (if \(g'\) is a duplication) or P3 (if \(g'\) is a speciation or a leaf), we conclude that \(\epsilon (g',s)=\textsf{True}\). For the second child, we have \(M_{T}(g'')\preceq s\), thus \(\textsf{F} _{T}_{V(Tg'')}\) is feasible for \((g'',s,X)\) and by P5, \(\delta ^\downarrow (g'',s)=\textsf{True}\). Finally, \(\delta ^*(g,s)=\textsf{True} = \delta (g,s)\).
(P2, \(\Rightarrow\)): Let \(\delta (g,s)=\textsf{Unknown}\). (Case P2.a) If \(s \in X\) then, from (1) g is internal and \(\delta ^*(g,s) = \textsf{Unknown}\). W.l.o.g., we may assume that \(\epsilon (g',s) \wedge \delta ^\downarrow (g'',s)=\textsf{Unknown}\). Since \(s \in X\), from (A1), we conclude that \(\delta ^\downarrow (g'',s)=\textsf{True}\). Next, \(\epsilon (g',s)=\textsf{Unknown}\) and by (A2) there is a weakly feasible mapping \(\textsf{F} _{T'}\) for \((g',s,X)\) and \(g'\) is an sduplication. Since \(s \in X\), similarly to the proof of P1, we can construct a feasible mapping for \((g',s,X)\) from weakly feasible \(\textsf{F} _{T'}\). A contradiction. (Case P2.b) Assume that \(s \notin X\). Then, from (2) g is internal and \(\delta ^*(g,s) \wedge \textsf{Unknown} = \textsf{Unknown}\). We have that \(\delta ^*(g,s) \in \{ \textsf{True}, \textsf{Unknown} \}\). (Case P2.b.1) Let \(\delta ^*(g,s)=\textsf{True}\). W.l.o.g. assume that \(\epsilon (g',s)=\delta ^\downarrow (g'',s) = \textsf{True}\). From P1 and \(s \notin X\), \(\delta (g',s)\) cannot be \(\textsf{True}\), thus, \(\sigma (g',s)=\textsf{True}\). From P3, there is a feasible \(\textsf{F} _{T'}\) for \((g',s,X)\) and \(\textsf{M} _{T'}(g')=s\). From P5, there is a feasible \(\textsf{F} _{T''}\) for \((g'',s,X)\). Let \(T=(T',T'')\), then \(\textsf{F} _{T'}\oplus \textsf{F} _{T''}\), is weakly feasible for (g, s, X). Note that there is no feasible mapping for (g, s, X), which follows from the fact that \(g'\) is an sspeciation, thus g is an sduplication and \(s \notin X\).
(Case P2.b.2) Let \(\delta ^*(g,s)=\textsf{Unknown}\). W.l.o.g. assume that \(\epsilon (g',s) \wedge \delta ^\downarrow (g'',s) = \textsf{Unknown}\). If \(\epsilon (g',s)=\textsf{True}\) then \(\delta ^\downarrow (g'',s)=\textsf{Unknown}\). Then, similarly to the previous case \(g'\) is an sspeciation and there is a feasible \(\textsf{F} _{T'}\) for \((g',s,X)\), while, from P6, there is a weakly feasible \(\textsf{F} _{T''}\) for \((g'',s,X)\). Then, similarly to the previous case, let \(T=(T',T'')\). Then, \(\textsf{F} _{T'}\oplus \textsf{F} _{T''}\) is weakly feasible mapping for (g, s, X). It remains to analyse the case when \(\epsilon (g',s)=\textsf{Unknown}\). By (A2) there is a weakly feasible \(\textsf{F} _{T'}\) for \((g',s,X)\) and \(g'\) is an sduplication. Here, \(\delta ^\downarrow (g'',s) \in \{\textsf{True},\textsf{Unknown} \}\), and depending on the value either, by P5 there is a feasible (\(\textsf{True}\)) or, by P6, weakly feasible (\(\textsf{Unknown}\)) \(\textsf{F} _{T''}\) for \((g'',s,X)\). We conclude that \(\textsf{F} _{T'}\oplus \textsf{F} _{T''}\) is weakly feasible mapping for (g, s, X). Since \(g'\) is an sduplication and \(s \notin X\), there is no feasible mapping for (g, s, X) in this case.
(P2, \(\Leftarrow\)). Let \(\textsf{F} _{T}\) be a weakly feasible mapping for (g, s, X) such that g is an sduplication in T. Since g is an sduplication, g is upper and \(s\notin X\). Thus, \(\delta (g,s)=\delta ^*(g,s) \wedge \textsf{Unknown}\) from (2). W.l.o.g. we may assume that \(\textsf{M} _T(g')=s\) (note that g is an sduplication). If \(g'\) is a speciation or a leaf, then \(\textsf{F} _Tg'\) is feasible for \((g',s,X)\) since no upper duplication is present in \(Tg'\). From P3, \(\sigma (g',s)=\epsilon (g',s)=\textsf{True}\). If \(g'\) is a duplication then \(g'\) is upper in \(Tg'\), thus \(\textsf{F} _Tg'\) is weakly feasible. From P2, \(\delta (g',s)=\textsf{Unknown}\). In all cases, \(\epsilon (g',s) \ge \textsf{Unknown}\). Similarly, for the second child of g, the mapping \(\textsf{F} _Tg''\) is either weakly feasible for \((g'',s,X)\) if \(g''\) is upper duplication outside X, or feasible otherwise. By P5 and P6, \(\delta ^\downarrow (g'',s) \ge \textsf{Unknown}\). Finally, \(\delta ^*(g,s) \ge \textsf{Unknown}\) and \(\delta (g,s)=\delta ^*(g,s) \wedge \textsf{Unknown} = \textsf{Unknown}\).
(P3, \(\Rightarrow\)). Let \(\sigma (g,s)=\textsf{True}\). Note that at least one of g and s is internal by the inductive assumption. Then, both nodes must be internal, from (7). W.l.o.g. we may assume that \(\delta ^\downarrow (g',s') \wedge \delta ^\downarrow (g'',s'')=\textsf{True}\). Thus, from P5, we have two feasible mappings \(\textsf{F} _{T'}\) for \((g',s',X)\) and \(\textsf{F} _{T''}\) for \((g'',s'',X)\). Let \(T=(T',T'')\), then g is an sspeciation and \(\textsf{F} _{T'} \oplus \textsf{F} _{T''}\) is feasible for (s, g, X).
(P3, \(\Leftarrow\)). Assume that is a feasible mapping \(\textsf{F} _T\) for (g, s, X) such that g is an sspeciation or an sleaf in T. If g is a leaf, the statement is obvious. Assume that g is an sspeciation, then s is internal and \(\sigma (g,s)\) follows from (7). W.l.o.g. we may assume that \(\textsf{M} _{T}(g') \preceq s'\) and \(\textsf{M} _{T}(g'') \preceq s''\). Thus, \(\textsf{F} _{T}_{V(Tg')}\) is feasible for \((g',s',X)\) since \(\textsf{F} _T(V(Tg'))\subseteq V(Ss')\). From P5, \(\delta ^\downarrow (g',s')=\textsf{True}\). Similarly, we obtain \(\delta ^\downarrow (g'',s'')=\textsf{True}\). Finally, \(\sigma (g,s) = \textsf{L}( \delta ^\downarrow (g',s') \wedge \delta ^\downarrow (g'',s'')) = \textsf{True}\).
(P4) It follows easily from the definition of \(\sigma\) and the operator \(\textsf{L}\).
(P5, \(\Rightarrow\)) Assume that \(\delta ^\downarrow (g,s)=\textsf{True}\). (Case P5.1) If s is a leaf, then g is internal by the inductive assumption. Then, by (4), \(\delta ^\downarrow (g,s)=\epsilon (g,s)=\delta (g,s) \vee \sigma (g,s) = \textsf{True}\). Note, that \(\sigma (g,s)=\textsf{False}\), otherwise both g and s must be leaves. Thus, \(\delta (g,s)=\textsf{True}\) and the feasible mapping for (g, s, X) exists by the already proven P1. (Case P5.2). If s is internal and \(s \notin X\), then, by (6), \(\delta (g,s) \vee \sigma (g,s) \vee \delta ^\downarrow (g,s') \vee \delta ^\downarrow (g,s'') = \textsf{True}\). If \(\delta (g,s)=\textsf{True}\), there is a feasible mapping for (g, s, X) from already proven P1 for g and s. Similarly, we have the mapping from P3 if \(\sigma (g,s)=\textsf{True}\). If \(\delta ^\downarrow (g,s')=\textsf{True}\), then there is a feasible mapping \(\textsf{F} _T\) for \((g,s',X)\) from the inductive assumption for P5. It should be clear that the mapping obtained from \(\textsf{F} _T\) by enlarging the codomain from \(V(Ss')\) to V(Ss) is feasible for (g, s, X). The remaining case when \(\delta ^\downarrow (g,s'')=\textsf{True}\) is analogous. (Case P5.3). If s is internal and \(s \in X\), then, by (5), \(\delta (g,s) \vee \sigma (g,s) \vee \textsf{M}\delta ^\downarrow (g,s') \vee \textsf{M}\delta ^\downarrow (g,s'') = \textsf{True}\). The proof is the same as above if \(\delta (g,s)=\textsf{True}\), \(\sigma (g,s)=\textsf{True}\), \(\delta ^\downarrow (g,s')=\textsf{True}\) or \(\delta ^\downarrow (g,s'')=\textsf{True}\). For the remaining case, let \(\delta ^\downarrow (g,s')=\textsf{Unknown}\), then there is a weakly feasible mapping \(\textsf{F} _T\) for \((g,s',X)\) from the inductive assumption for P6. By placing all upper duplications from \(\textsf{F} _T\) at episode s and enlarging the codomain from \(V(Ss')\) to V(Ss), we construct a feasible mapping for (g, s, X). The remaining case when \(\delta ^\downarrow (g,s'')=\textsf{Unknown}\) is analogous.
(P5, \(\Leftarrow\)) Assume there is a feasible mapping for (g, s, X). (Case P5.1) If s is a leaf, then g is internal by the inductive assumption. Thus, g is an upper sduplication in T. Thus, \(s \in X\) and by already proven P1, \(\delta (g,s)=\textsf{True}\), \(\epsilon (g,s)=\textsf{True}\) and \(\delta ^\downarrow (g,s)=\textsf{True}\). (Case P5.2) Assume s is an internal node. If g is an sduplication, then similarly to the above case, from P1, \(s \in X\), and \(\delta (g,s)=\textsf{True}\) and \(\delta ^\downarrow (g,s)=\textsf{True}\) using (5). Analogously, we have the same conclusion, if g is an sspeciation (here, s does not have to be in X). For the remaining cases, there is \(v \prec s\), such that g is either an vduplication, or an vspeciation. W.l.o.g., we may assume that \(v \preceq s' \prec s\). (Case P5.2.a) If g is an vspeciation, then \(\textsf{F} _T(g)=v\) and there is no upper duplication in T. Since Tv is a subtree of \(Ts'\), \(\textsf{F} _T\) with a codomain \(V(Ts')\) is feasible for \((g,s',X)\). From, the inductive assumption for P5, \(\delta ^\downarrow (g,s')=\textsf{True}\) and also \(\textsf{M}\delta ^\downarrow (g,s') = \textsf{True}\) (if \(s \in X\)). In all cases, \(\delta ^\downarrow (g,s)=\textsf{True}\). (Case P5.2.b) If g is an vduplication, then, from the feasibility of \(\textsf{F} _T\), there is at least one node \(w \in X\) on the path from s to v (there must be a duplication episode for g). Let w be the lowest (i.e., closest to v) node with the property. Now, we have two cases. If \(w=s\), then \(s \in X\), and there is no candidate in X for g below s, therefore, there is no the feasible mapping for \((g,s',X)\). However, there is a weakly feasible \(\textsf{F} '_T\) for \((g,s',X)\) infered from \(\textsf{F} _T\) by setting \(\textsf{F} '_T(d):=s'\) all upper duplications d in T having \(\textsf{F} _T(d)=s\), and \(\textsf{F} '_T(u):=\textsf{F} _T(u)\) for the remaining nodes u from T. Clearly, \(\textsf{F} '_T(V(T)) \subseteq V(Ss')\). Hence, from the inductive assumption P6, \(\delta ^\downarrow (g',s)=\textsf{Unknown}\). Next, \(\textsf{M}\delta ^\downarrow (g',s) = \textsf{True}\) and \(\delta ^\downarrow (g,s)=\textsf{True}\) using (5). It remains to analyse the case when, \(w \preceq s' \prec s\). Since, \(\textsf{F} _T\) is \(\preceq\)monotonic, \(\textsf{F} _T(V(T)) \subseteq V(Sw) \subseteq V(Ss')\). Then, \(\textsf{F} _T\) with the codomain shrinked to \(V(Ss')\) is feasible for \((g,s',X)\). From the inductive assumption for P5, \(\delta ^\downarrow (g,s')=\textsf{True}\). Here, it does not matter whether \(s \in X\), in both cases we get \(\delta ^\downarrow (g,s)=\textsf{True}\) from (5) or (6).
(P6, \(\Rightarrow\)) Assume, \(\delta _X^\downarrow (g,s)\) is \(\textsf{Unknown}\). Here, \(\textsf{Unknown}\) is obtained only from (6), i.e., when \(s \notin X\) and s is internal. We have \(\delta (g,s) \vee \sigma (g,s) \vee \delta ^\downarrow (g,s') \vee \delta ^\downarrow (g,s'') = \textsf{Unknown}\). By P4, \(\sigma (g,s)=\textsf{False}\). If \(\delta (g,s)=\textsf{Unknown}\), then, we have a weakly feasible mapping from already proven P2. If \(\delta ^\downarrow (g,s')=\textsf{Unknown}\), then, we have a weakly feasible mapping from the inductive assumption for P6. The same holds for \(s''\).
(P6, \(\Leftarrow\)) Assume there is a weakly feasible mapping \(\textsf{F} _T\) for (g, s, X). (Case P6.1) If s is a leaf, then g is internal by the inductive assumption. Thus, g is an upper sduplication in T. Thus, \(s \notin X\) and by already proven P2, \(\delta (g,s)=\textsf{Unknown}\). Since, \(\sigma (g,s)=\textsf{False}\), \(\epsilon (g,s)=\textsf{Unknown}\) and \(\delta ^\downarrow (g,s)=\textsf{Unknown}\) from (4). (Case P6.2) Assume s is internal. Note that g cannot be an vspeciation for any v, otherwise there is no upper duplication in T and \(\textsf{F} _T\) cannot be weakly feasible. If g is an sduplication in T, then similarly to the previous case, from P2, \(s \notin X\), and \(\delta (g,s)=\textsf{Unknown}\). Note that \(\delta ^\downarrow (g,s')=\delta ^\downarrow (g,s'')=\textsf{False}\) (since \(\textsf{M} _T(g)=s\)). We conclude that \(\delta ^\downarrow (g,s)=\textsf{Unknown}\) from (6). For the remaining case, there is \(v \prec s\), such that g is an vduplication. W.l.o.g., we may assume that \(v \preceq s' \prec s\). From the weak feasibility of \(\textsf{F} _T\), no node from the path between s to v (inclusively) is in X. Similarly, to P5.2, we construct weakly feasible \(\textsf{F} '_T\) for \((g,s',X)\) from \(\textsf{F} _T\) by setting \(\textsf{F} '_T(d):=s'\) all upper duplications d in T having \(\textsf{F} _T(d)=s\), and \(\textsf{F} '_T(u):=\textsf{F} _T(u)\) for the remaining nodes u from T. By the inductive assumption for P6, \(\delta ^\downarrow (g',s)=\textsf{Unknown}\). Next, if \(v=s'\), then \(\delta (g,s')=\textsf{Unknown}\) (since \(x \notin X\)) by P2, otherwise \(\delta (g,s')=\textsf{False}\). Also, \(\delta ^\downarrow (g,s'')=\textsf{False}\). This yields \(\delta (g,s)=\textsf{Unknown}\) using (6). \(\quad\quad\quad \square\)
To solve Episode Feasibility, we have to apply \(\delta ^\downarrow\) on the roots of the input trees.
Theorem 2
(Correctness) Given a partial gene tree G over a species tree S and \(X \subseteq V(S)\). G is Xfeasible if and only if \(\delta _X^\downarrow (\textsf{root} (G),\textsf{root} (S))\) is \(\textsf{True}\).
Proof
The proof follows immediately from P5 of Lemma 1: \(\delta _X^\downarrow (\textsf{root} (G),\textsf{root} (S))\) is \(\textsf{True}\) if and only if there is a feasible \(\textsf{F} _T\) for \((\textsf{root} (G),\textsf{root} (S),X))\) such that T extends G and \(\textsf{F} _T(\textsf{Dup} _T) \subseteq X\). \(\quad\quad \square\)
Theorem 3
(Complexity) Given a partial gene tree G over a species tree S and \(X \subseteq V(S)\). The time and space complexity of solving Episode Feasibility by the dynamic programming algorithm is O(V(G)V(S)).
Proof
We have three arrays \(\delta _X\), \(\delta ^\downarrow _X\) and \(\sigma _X\) (note that \(\epsilon\) and \(\delta ^*\) can be directly inserted in their calls), each of size O(V(G)V(S)) and every cell of an array can be computed in O(1) time. \(\quad\quad\quad\quad\quad \square\)
An example of DP execution with a feasible solution is depicted in Fig. 2.
Solving MetaEC for a single partial gene tree
Here we describe the main algorithm to solve MetaEC for instances with a single gene tree. First, we characterize an important property of episodes.
Lemma 4
(Fixed Episodes) Given a partial gene tree G over a species tree S. Assume that there are nodes g in G and s in S such that

if s is the root of S, then at least one proper subtree of G contains species (leaflabels) from both children of s.

otherwise, let p be the parent of s, then Gg is a gene tree, g is a pspeciation and a child of g is an sduplication.
Then, for any \(G^*\) that extends G, s is an episode in every valid mapping between \(G^*\) and S.
Proof
For the first case, all nodes above the root of the subtree are gene duplications mapped to the root of S in any \(G^*\) that extends G. Therefore, the root of S is an episode in all valid mappings. In the second case, the duplication child cannot be raised, therefore, its mapping is fixed. \(\square\)
The nodes satisfying the above conditions we call fixed episodes (for G and S). For example, for trees from Fig. 1, there are two fixed episodes: the root of S and the leaf b, where the duplications with fixed mappings are depicted using white marks in the exemplary gene tree G. The set of all fixed episodes can be computed in linear time and space by bottomup traversal of the partial gene tree G and by using LCAqueries in the species tree S as follows. For each node g from V(G), the algorithm computes a tuple (u, s, d), where

\(u \in \{\textsf{True},\textsf{False} \}\) is \(\textsf{True}\) if and only if there is a leaf with unknown label reachable from g,

\(s \in V(S) \cup \{\textsf{None}\}\) is the least common ancestor of all non\(\bot\) labels reachable from g in S and \(\textsf{None}\) if only \(\bot\)’s are visible from g,

and \(d \in \{\textsf{True},\textsf{False} \}\) is \(\textsf{True}\) if and only if \(u=\textsf{False}\) and g is a duplication node in a gene tree Gg.
Then, for each g and its tuple (u, s, d), and for each child of g with a tuple \((u',s',d')\):

if \(s=s'=\textsf{root} (S)\), then s (the root of S) is a fixed episode,

if \(u=d=\textsf{False}\), \(d'=\textsf{True}\) and the parent of \(s'\) is s, then \(s'\) is a fixed episode.
We omit correctness and complexity proofs for brevity. Note that the number of fixed episodes is the lower bound of \(\textsf{EC} (G,S)\).
Algorithm 1 takes as input a partial gene tree G over a species tree S and outputs \(\textsf{EC} (G,S)\). It first computes the set of fixed episodes F (see Lemma 4). The algorithm then starts with an initial maximal episode number b equal to the number of nodes in S. In each iteration of a while loop, the algorithm checks if there is a set C of size \(bF1\) from the vertices of S that are not in F, such that the partial gene tree G is \(C \cup F\)feasible using the dynamic programming algorithm. This step requires \(\left( {\begin{array}{c}V(S)F\\ bF1\end{array}}\right)\) calls of DP in the worst case. If such a set C exists, the algorithm computes \(\textsf{EC} (G^*,S)\) by the linear time algorithm from [20], where \(G^*\) is the gene tree obtained by backtracking from the corresponding call of DP and updates b with the result. Note that b is not assigned the value of \(C \cup F\), since the minimal set of episodes for \(G^*\) and S is a subset of \(C \cup F\), and it is often significantly smaller than \(C \cup F\) in early steps of iteration. Updating b with \(\textsf{EC} (G^*,S)\) guarantees the minimal number of episodes, where some elements of C may be unused. This is an important optimization step. If such a set C does not exist, the algorithm terminates and returns the current value of b.
The correctness of the algorithm follows from the fact that if there is no set X of size \(b1\) such that G is Xfeasible, then there is no set of any size smaller than b that satisfies the property. Since b represents the number of episodes from some valid mapping, it is also minimal in such a case. Therefore, when the algorithm terminates, \(b=\textsf{EC} (G,S)\), and the algorithm returns the correct value.
The algorithm’s worstcase time complexity is \(\sum _{k=f}^{nf}{\left( {\begin{array}{c}nf\\ k\end{array}}\right) nm} = O(nm2^n)\), where f is the size of the set of fixed episodes (\(f=F\)), n denotes the number of vertices in S, and m denotes the number of vertices in G. Despite the exponential time complexity, in our experiments on both simulated and empirical data, we were able to compute exact solutions after only a few executions of the main loop.
Distributions of genespecies mappings
To evaluate the accuracy of genespecies mappings, we propose a method that enhances the DP algorithm by incorporating formulas for inferring the number of genespecies mappings present in all feasible reconstructed mappings. These counts can be collectively integrated to determine, for any leaf with unknown label, the precise frequency of its mappings to each species leaf within these feasible mappings. An alternative method for approximating these frequencies was suggested in [31] using uniform sampling. However, in this work, we introduce an exact algorithm for this purpose.
For a fixed species tree S, and a node g from a partial gene tree G, we call a mapping \(f :L(Gg) \times L(S) \rightarrow \{0,1,2,\dots \}\) a counter on L(Gg) if for every leaf l with unknown label, the sum \(\sum _{s} f(l,s)\), denoted \(\#f\), does not depend on l. Counters will be used to count how many times a given gene \(\bot\)leaf is assigned to a species leaf in all feasible mappings. In such a case, \(\# f\) is the number of all such feasible mappings. The counter fixed to a gene tree \(\bot\)leaf l and represented by the function \(f(l,\cdot )\), is referred to as the genespecies distribution (of l). Subsequently, we often examine normalized distributions, wherein each value is divided by \(\# f\) (then \(\sum _s f(l, s)=1\)). For convenience counters also include other leaves, but their counts will be set to 0. We have the following basic counters:

\(\emptyset _A\) is the zero counter on A, i.e., the counter with \(\# \emptyset = 0\),

for s in L(S), \(\textsf{B} _{l,s}\) is a counter on \(\{l\}\) such that \(\textsf{B} _{l,s}(l,s)=1\), and \(\textsf{B} _{l,s}(l,s')=0\) for all \(s' \ne s\).
Let \(\oplus\) and \(\otimes\) be commutative operators, where \(\otimes\) has higher precedence that \(\oplus\), satisfying the following properties.

If f and g are counters on disjoint sets A and B, respectively, then \(f \otimes g\) is the counter on \(A \cup B\), such that for every \(l \in A\), \((f \otimes g)(l,s)=f(l,s) \cdot \#g\).

If f and g are counters on A, then \(f \oplus g\) is the counter on A, such that, for every \(l \in A\), \((f \oplus g)(l,s)=f(l,s) + g(l,s)\).
The empty counter is used when a part of a gene tree has no leaves with unknown labels, while the counter \(B_{l,s}\) represents the situation when leaf l is assigned to s. Additionally, we use a special counter \(\textsf{E}\) that represents situations where DP formulas return \(\mathbb {F}\). The counter is evaluated as follows: \(\textsf{E} \oplus f = f\) and \(\textsf{E} \otimes f = \textsf{E}\) for any counter f including \(\textsf{E}\).
By using the notation introduced in the dynamic programming algorithm, we define counters for DP \(\delta _{\#\mathbb {T}}(g,s)\), \(\delta _{\#\mathbb {U}}(g,s)\), \(\delta ^\downarrow _{\#\mathbb {T}}(g,s)\), \(\delta ^\downarrow _{\#\mathbb {U}}(g,s)\), \(\epsilon _{\#\mathbb {T}}(g,s)\), \(\epsilon _{\#\mathbb {U}}(g,s)\) and \(\sigma _{\#}(g,s)\) to count distributions of genespecies mappings. E.g., \(\delta _{\#\mathbb {T}}(g,s)\) is a counter for the leaves from L(Gg) that counts mappings for these leaves for the cases when \(\delta (g,s)\) is \(\textsf{True}\), and so on.

In (1), \(\delta _{\#\mathbb {T}}(g,s) = e_{\mathbb {T},\mathbb {T}} \oplus r\) where, \(e_{i,j}=\epsilon _{\#i}(g', s) \otimes \epsilon _{\#j}(g'', s)\) and if \(\textsf{sib} (v)\) is the sibling of a node v, \(r=\bigoplus _{i,j \in \{\mathbb {T},\mathbb {U}\}, v \in \{g',g''\}, w \in \{s',s''\}} \epsilon _{\#i}(v, s) \otimes \delta ^\downarrow _{\#j}(\textsf{sib} (v),w)\).

In (1–2) \(\delta _{\#\mathbb {U}}(g,s) = e_{\mathbb {T},\mathbb {T}} \oplus e_{\mathbb {U},\mathbb {T}} \oplus e_{\mathbb {T},\mathbb {U}} \oplus e_{\mathbb {U},\mathbb {U}} \oplus r\), where, r and \(e_{i,j}\) are defined above.

In (10), \(\epsilon _{\#\mathbb {T}}(g,s) = \sigma _{\#}(g,s) \oplus \delta _{\#\mathbb {T}}(g,s)\), and, if \(\delta (g,s)=\textsf{Unknown}\), then \(\epsilon _{\#\mathbb {U}}(g,s) = \delta _{\#\mathbb {U}}(g,s)\).

In (4), \(\delta ^\downarrow _{\#\mathbb {U}}(g,s)=\epsilon _{\#\mathbb {U}}(g,s)\) and \(\delta ^\downarrow _{\#\mathbb {T}}(g,s)=\epsilon _{\#\mathbb {T}}(g,s)\).

In (5), \(\delta ^\downarrow _{\#\mathbb {T}}(g,s)= \epsilon _{\#\mathbb {T}}(g,s) \oplus \delta ^\downarrow _{\#\mathbb {T}}(g,s') \oplus \delta ^\downarrow _{\#\mathbb {U}}(g,s') \oplus \delta ^\downarrow _{\#\mathbb {T}}(g,s'') \oplus \delta ^\downarrow _{\#\mathbb {U}}(g,s'')\) and \(\delta ^\downarrow _{\#\mathbb {U}}(g,s)=\epsilon _{\#\mathbb {U}}(g,s)\).

In (6), \(\delta ^\downarrow _{\#\mathbb {T}}(g,s)= \epsilon _{\#\mathbb {T}}(g,s) \oplus \delta ^\downarrow _{\#\mathbb {T}}(g,s') \oplus \delta ^\downarrow _{\#\mathbb {T}}(g,s'')\) and \(\delta ^\downarrow _{\#\mathbb {U}}(g,s)= \epsilon _{\#\mathbb {U}}(g,s) \oplus \delta ^\downarrow _{\#\mathbb {U}}(g,s') \oplus \delta ^\downarrow _{\#\mathbb {U}}(g,s'')\).

In (7), \(\sigma _\#(g,s) = \delta ^\downarrow _{\#\mathbb {T}}(g',s') \otimes \delta ^\downarrow _{\#\mathbb {T}}(g'',s'') \oplus \delta ^\downarrow _{\#\mathbb {T}}(g',s'') \otimes \delta ^\downarrow _{\#\mathbb {T}}(g'',s')\).

In (8), \(\sigma _\#(g,s) = \textsf{B} _{g,s}\) if \(g = \bot\), otherwise \(\sigma _\#(g,s) = \emptyset _{\{g\}}\).

In all remaining uncovered cases the counters are equal to \(\textsf{E}\).
The following lemma states the crucial property of DP counters.
Lemma 5
(Correctness of counters for DP) Given a partial gene tree G and a species tree S and \(X\subseteq V(S)\). G is Xfeasible if and only if the counter \(\delta ^\downarrow _{\#\mathbb {T}}(\textsf{root} (G),\textsf{root} (S))\) is not \(\textsf{E}\) and for every leaf l of G with unknown label, \(\delta ^\downarrow _{\#\mathbb {T}}(\textsf{root} (G), \textsf{root} (S))(l,s)\) is the number of all gene trees \(G^*\) extending G such that \(\textsf{F} _{G^*}\) is a valid mapping satisfying \(\textsf{F} _{G^*}(\textsf{Dup} _{G^*}) \subseteq X\) and \(\textsf{F} _{G^*}(l)=s\).
The proof of the above lemma follows by induction, similar to the proof of correctness of DP. We omit technical details.
An example of genespecies distributions is depicted in Fig. 3. DP counters with verification algorithm are implemented in the software package metaEC.
Extensions
To identify the optimal solution within the main loop, enumerating all possible combinations of size \(bf1\) from the set of episode candidates \(V(S)\setminus F\) may be timeconsuming for larger instances. To address this issue, we propose a heuristic approach that randomly samples combinations of size \(bf1\) if \(\left( {\begin{array}{c}nf\\ bf1\end{array}}\right)\) is large (e.g., \(>1000\)) and adds a stopping condition based on the number of dynamic programming (DP) calls without improvement (e.g., after 100 calls). This approach not only speeds up the algorithm but also provides additional information on whether the returned value is exact or an upper bound obtained by switching to a heuristic mode. See Sect. 4 for more details.
Furthermore, based on our experiments, we have observed that the solution is often close to the set of fixed episodes. To leverage this observation, we propose a bottomup algorithm that explores candidate sets starting from sizes 0, 1, 2, and so on until a feasible solution is found. In this case, the internal search has a time complexity of \(O(\left( {\begin{array}{c}nf\\ i\end{array}}\right) )\), starting from \(i=0\). This algorithm can be combined with the heuristic variant described earlier to improve its effectiveness. However, the experimental evaluation did not show significant improvement compared to the topdown method in Algorithm 1.
MetaEC in the general case
Here we show that MetaEC in a general case can be solved using a single partial gene tree under an additional assumption. Given a collection of partial gene trees \(G_1, G_2, \dots , G_k\) over a species tree S, let \(\omega\) be a new species, called outgroup, not present in S. We first add the outgroup to every input tree. Let \(S^\omega\) be species tree \((S,\omega )\), \(G_1^\omega =(G_1,\omega )\) and \(G^\omega _i=((G_i,\omega ),G^\omega _{i1})\), for \(i>1\). Then, by \(\omega\)MetaEC we define the problem MetaEC with a single partial gene tree, where the extension of a partial labeling cannot introduce \(\omega\), i.e., if \(\Lambda _{G_1}(v)=\bot\) then \(\Lambda _{G^*_1}(v) \ne \omega\). See Fig. 4 for illustration.
We have the following property.
Lemma 6
Given a collection of at least two partial gene trees \(G_1, G_2, \dots , G_k\) over a species tree S such that \(\omega \notin L(S)\). \(X \subseteq V(S)\) is the set of episodes that yields the solution of MetaEC for \(G_1, G_2, \dots , G_k\) and S if and only if \(X \cup \{\textsf{root} (S^\omega )\}\) is the set of episodes that yields the solution to the instance \(G^\omega _k\) and \(S^\omega\) of \(\omega\)MetaEC.
Proof
(\(\Rightarrow\)) Assume that \(G^*_i\) extends \(G_i\), for each i. Then, \((G^*_i)^\omega\) extends \(G_i^\omega\). Since the extension of \(G_i\) introduces only nodes from L(S), we have the same property with \((G^*_i)^\omega\). Now, every parent of a leaf labeled \(\omega\) in \((G^*_k)^\omega\) is speciation mapped to the root of \(S^\omega\), since, for some i, its sibling is a root of \(G^*_i\) and \(\omega \notin \mathcal {L} (G^*_i)\). Thus, if \(i>1\), the root of \((G^*_i)^\omega\) is a duplication mapped to the root of \(S^\omega\). In summary, all duplications from gene trees \(G^*_i\) in \((G^*_k)^\omega\) are mapped below the root of \(S^\omega\), and they are separated by speciation nodes from the duplications mapped to the root of \(S^\omega\) as indicated in Fig. 4. Now, we define a valid mapping \(F^\omega _k = \textsf{F} _{(G^*_k)^\omega } :V((G^*_k)^\omega ) \rightarrow V(S^\omega )\). For each i, \(\textsf{F} ^\omega _k\) on the set of nodes of \(G_i\) equals the corresponding valid mapping between \(V(G^*_i)\) and V(S) that yields the solution to MetaEC, while for the remaining nodes v we set \(\textsf{F} _k^\omega (v) = \textsf{root} (S^\omega )\). It should be clear that \(F^\omega _k\) is a valid mapping. Now, it is not difficult to see that the set of duplication episodes in \((G^*_k)^\omega\) is \(X \cup \textsf{root} (S^\omega )\). If there is a solution to \(\omega\)MetaEC with a lower number of episodes than \(X+1\), say obtained by \(X' \cup \{\textsf{root} (S^\omega )\}\) with \(X'<X\), then the construction can be reversed to obtain valid mappings and the corresponding duplication episodes \(X'\) for the initial instance of MetaEC. However, this is a contradiction with the assumption that X is the solution to the initial instance of MetaEC.
(\(\Leftarrow\)) The proof of the second direction is analogous since the transformation between collections of partial gene trees and the partial gene tree \(G^\omega _k\) is reversible. We omit easy details. \(\quad\quad\quad\quad\quad \square\)
Note that the algorithms provided in the previous sections can be easily modified to solve \(\omega\)MetaEC, by replacing case (8) with:
Then, DP will exclude extensions of \(\bot\) by \(\omega\).
Experiments
In this Section we present two computational studies based on simulated and empirical data.
Simulated dataset
The species trees having none, one, or two whole genome duplication events were taken directly from [25]. Then, we estimated gene trees via tree inference software from simulated sequences and modified them to represent the uncertainty commonly associated with metagenomic data. Finally, our algorithm was evaluated on six datasets consisting of estimated gene trees.
To begin, we describe how the species and gene trees were generated. Then, we explain the modifications made to the gene trees to represent the uncertainty of metagenomic data. Finally, we show how the results obtained with our algorithms allow us to infer genomewide duplication events and gene and species distributions.
A species tree
First, we briefly summarize the simulation procedure from [25]. The simulated species trees were generated by SimPhy [33] with parameter settings used in a simulated study [34] that was based on an empirical dataset of 16 Fungi species [35]. The species tree S of 20 taxa was generated by SimPhy with the speciation rate parameter equal to \(1.8 \times 10^{9}\) and the tree height parameter set to \(1.8 \times 10^9\).
To simulate a whole genome duplication (WGD), a node v in the species tree S was chosen as the location of the event. Subsequently, a modified species tree, denoted as \(S'\), was constructed by substituting a subtree Sv with a duplicated version of itself. This duplication involved creating a new root connected to the original root of Sv and the root of its copy. The WGD variants used in the simulations are illustrated in Fig. 5, where \(S_1\) represents a single recent WGD \(\Gamma\), \(S_2\) represents a single ancient WGD \(\Theta\), \(S_3\) represents two WGDs \(\Lambda\) and \(\Delta\), with \(\Delta\) occurring after \(\Lambda\), \(S_4\) represents two close WGDs \(\Pi\) and \(\Pi '\) at the same branch, and \(S_5\) represents two recent independent WGDs \(\Psi\) and \(\Phi\). Let us refer to the original tree S with no WGD events as \(S_0\). Note that these simulation methods do not incorporate fractionation, that is, the loss of a gene copy to eliminate the redundancy [36].
Gene trees
For every \(S_i\), \(i \in \{0, 1, \dots , 5\}\), one hundred true gene trees were generated using SimPhy. The duplication and loss rate parameter was set to \(2^{10}\) events per generation per lineage. To minimize the effect of incomplete lineage sorting, the population size parameter was set to 10. All other parameters were taken from [34].
Next, we describe our pipeline to infer estimated trees from true trees. For every true tree G, we first simulated a multiple sequence alignment (MSA).
For the MSA simulation, we used INDELible [37] and parameters from [34], with one difference, we used a constant sequence length of 1000. Then, we inferred an unrooted maximumlikelihood tree from each MSA using FastTree [38] with the GTR model. Finally, we performed midpointplateau rooting of each unrooted gene tree using URec [39]. The rooting was inferred in a way that minimizes the duplicationloss cost between the gene tree and the species tree S.
In summary, we obtained six datasets of estimated rooted gene trees denoted \(\mathcal G_i=\{G_{i,1},G_{i_2},\dots ,G_{i,{100}}\}\), where each dataset comprises 100 gene trees generated using the same species tree S but with a different WGD scenario \(S_i\). In our evaluation study, for each dataset \(\mathcal G_i\), we generated a set of partial gene trees \(\tilde{\mathcal G}^{k,p}_i\) by randomly removing each leaf label from every gene tree \(G_{i,j}\) in \(\mathcal G_i\) with probability p. We considered values of p from 0.0 to 0.6 in increments of 0.1 and generated 100 instances of \(\tilde{\mathcal G}^{k,p}_i\) for each value of p and \(k=1,2,\dots ,K\), with \(K=100\). This resulted in a total of 4200 instances \((\tilde{\mathcal G}^{k,p}_i,S)\) of MetaEC, where S is the species tree used to generate \(\mathcal G_i\), and \(\tilde{\mathcal G}^{k,p}_i\) consists of 100 partial gene trees. Note that the instances for \(p=0\), correspond to no removal of leaf labels, i.e., \(\tilde{\mathcal G}^{k,0.0}_i=\mathcal G_i\), for each i and k.
WGD detection
The results of our algorithm on the estimated partial gene trees are depicted in Fig. 5, where we summarize the episode sizes in the form of histograms. Also, additional data is provided in Figs. 6 and 7. The evaluation took about 24 h of a computing server with 80 cores. In general, we observed that the runtime and the number of DP calls grow linearly with the value of parameter p on average as indicated in the bottom diagram of Fig. 8. We used the heuristic variant of the algorithm, where random sampling was applied if the number of combinations exceeded 100 trees, and with the stopping criterion equal to 50. Out of 4200 instances, 3037 were completed with the exact solution (see Fig. 6). The resulting costs without exact guarantee, were more often obtained for larger values of p. Additionally, we observed that the lower bound given by the number of fixed episodes was a tight approximation of the inferred cost (see Fig. 8). As p increased, the number of fixed episodes decreased. Note that for \(p=1\) (not included in our analysis), the solution to MetaEC is 1. In such a case the leaf labelings are constant functions and all internal nodes of each gene tree are duplications. Then there is just one episode, which can be placed at the root of the species tree.
Now, we briefly summarize the outcome of WGD detection. Initially in the preliminary version of the article [40], we quantified support by counting the number of single gene duplications aggregated within a particular episode. Here, to reflect the contributions of single gene families while mitigating the influence of large gene families with potentially numerous duplications, we count the number of gene trees that contain duplications clustered within the episode. Note that in contrast to [40] where the true trees were used to detect WGD events, here we perform the detection on estimated gene trees.
In the context of the scenario without WGDs (\(S_0\)), the episodes can be treated as background noise consisting of single duplication events. Consequently, when subtracting the contribution of background duplication episodes from the duplication clustering results on simulated datasets containing WGDs, we can emphasize the simulated WGD events as significant occurrences. Formally, our analysis is conducted using the following formulas.
Given a dataset \((\tilde{\mathcal G}^{k,p}_i,S)\) and a node s in S, the psupport of s (w.r.t. k and i), denoted \(\mu _{p,s,i,k}\) is the number of gene trees obtained by applying the DP algorithm on trees from the dataset \({\mathcal G}^{k,p}_i\), that have at least one duplication in the episode s. Then, the normalized average psupport of a species node s in scenario \(S_i\) is defined as
Recall that \(K=100\) is the number of estimated gene tree datasets for fixed p and i. In other words, the normalized average support accounts for the difference between the average support values of \(S_i\) and \(S_0\). For example, given that every gene tree dataset consists of 100 estimated gene trees, the normalized average psupport of s in the scenario \(S_i\) close to 100, denotes high support for the duplication episode at s (in \(S_i\)). Note also that the normalized psupport can be negative due to normalization.
The results for \(S_1\) (as depicted in Fig. 5) reveal that WGD \(\Gamma\) is uniquely wellsupported when \(p \le 0.4\), with support values exceeding 20. It is noteworthy that the outcomes obtained for \(S_1\) outperform the analysis conducted in [40], which utilized true trees, in terms of accurately identifying the simulated WGD event. In summary, the WGD detection outcome depends on the location of a WGD and the value of p. In general, when p increases, the support of WGDs decreases. Recent WGDs like \(\Gamma\), \(\Delta\), and \(\Phi\) are well supported up to \(p=.3\), while more ancient \(\Psi\) up to \(p=.2\). For ancient WGDs \(\Theta\), \(\Lambda\), \(\Pi\), \(\Pi '\), the normalized average psupport is low and, for most cases, indistinguishable from nonWGD nodes. However, the root episode is generally well supported in scenarios having ancient WGDs, which suggest such an event close to the root location.
Genespecies distributions
We implemented the algorithm outlined in Sect. Distributions of genespecies mappings to derive counters with genespecies distributions for \(\bot\)leaves across all gene tree datasets, i.e., for all six species trees \(S_0\)\(S_5\), all six different positive values of p (i.e., we inferred \(36 \times 100\) counter collections for each gene tree set). The results frequently yielded remarkably high numbers of feasible mappings, often exceeding \(10^{1000}\). As a remedy, we employed normalization by dividing the values of the genespecies distributions by \(\#f\), where f was the relevant counter.
We first analyzed the domain of distributions with positive values. Let \(d :L(S) \rightarrow [0,1]\) be a genespecies distribution. The span of d is the set of all species leaves with a positive value. Furthermore, if the span of this distribution is the set of leaves originating from a subtree rooted at s within the species tree S, we categorize this distribution as a subtreespanning distribution. A summary of subtreespanning and nonsubtreespanning distributions can be found in Table 2. It is noteworthy that spanning distributions occur with high frequency. Furthermore, they often encompass significant portions of the species tree, and the property gets stronger with the increase of p. This, however, is not the desired characteristic, as the most soughtafter distributions are the onepoint distributions where a single species leaf has the maximum value of 1.
Subsequently, we conducted an analysis to determine the extent to which these distributions deviate from uniformity. Detailed insights can be found in Fig. 9. In this figure, we present a computed coefficient of variation (CV) for each distribution, which is the standard deviation normalized by the mean. Lower CV values signify a closer resemblance to a uniform distribution. The majority of histograms display a tall blue bar indicating low CV values. This suggests that the distributions tended to exhibit a significant degree of nearuniformity in nearly all simulated datasets.
In summary, our study reveals that the species presented in the reconstructed gene leaf mappings extend across a substantial portion of the species tree. Additionally, the frequency distribution across all feasible mappings typically exhibits an almost uniform shape. As a result, identifying the true genespecies mapping signal within the leaf mapping proves challenging with the current approach. This challenge is, in part, attributed to the combinatorial explosion in the number of feasible mappings, consequently leading to flattened distributions of genespecies mappings based on leaf counts.
Empirical evaluation
To ensure that our algorithm was properly tested, we required a dataset that would capture the characteristics of the metagenomic data as closely as possible, while allowing us to assess the quality and accuracy of the results obtained. For this reason, we decided to prepare a dataset consisting of gene trees for species identified during metagenomic analysis. To simulate unknown genespecies assignments, we artificially removed some of the gene labels from the gene trees and retained information about their taxonomic origin for further analysis of the results. Another important issue was the presence of a previously described wholegenome duplication event that occurred in the evolutionary tree of selected species. Given the above requirements, we decided to use proteomes belonging to yeast species identified during metagenomic analysis of kefir [41]. Note that a direct comparison of our results with alternative methods is not possible due to the absence of any existing approach for the simultaneous inference of duplications and genespecies mappings.
Data preparation
The eight selected species are: Kazachstania Africana, Kazachstania naganishii, Naumovozyma dairenensis, Tetrapisispora blatte, Tetrapisispora phaffi, Torulaspora delbrueckii, Zygosaccharomyces rouxii and Saccharomyces cerevisiae. A species tree containing the listed species, consistent with the NCBI taxonomy and many papers on yeast evolution, is shown in Fig. 10. It also shows the location of the wholegenome duplication event confirmed by previous studies [42, 43].
The proteomes used to infer gene trees were sourced from the UniProt database [44]. Protein families were created by dividing the proteins into groups using the mcl program [45] with parameters \(I=2\) and \(I=5\). However, since the differences between the obtained sets were minimal, we used the set obtained for \(I=2\) in subsequent steps. The protein sequences in each group were aligned with the MUSCLE algorithm [46], and unrooted gene trees were inferred using the phyml program [47] with the default parameter setting. Rooting of the gene trees was performed by URec program [39] using the minimal duplicationloss cost as the rooting criterion.
We removed trees containing fewer than 3 leaves or 3 species, as well as trees with edges of length 0 from the final set of trees. This resulted in 3430 rooted gene trees. Similar to the first experiment, we created 10 datasets for each \(p \in \{0.1, 0.2, \dots , 0.6\}\) by randomly removing each leaf label with the probability p. This resulted in 60 datasets plus the original dataset representing \(p=0\).
Results
Figure 10 depicts histograms showing the results obtained for the described dataset. The evaluation was performed on the same computing server as before and took approximately one hour. For this evaluation, we set the sampling threshold and stopping criterion to 50, which yielded exact solutions for all cases. The number of fixed episodes was consistent across datasets with \(p<0.6\), at 12 (note that the number of leaves in a tree was 15). For \(p=0.6\), the number of fixed episodes fell within the range of 9 to 12. The number of DP calls ranged from 2 to 3 for \(p<0.6\) and between 2 and 11 for \(p=0.6\).
The results obtained by the algorithm for the yeast dataset are consistent with our knowledge of the wholegenome duplication localization. For the dataset with all leaf labels present and for \(p=0.1\), we have the highest support for the WGD event. The number of supporting single duplications decreases gradually for successive \(p's\). For values of \(p \ge 0.5\), the correct WGD localization is still supported by a significantly large number of single duplications. It is worth noting that even for the \(p=0.6\), the right location is supported by three times as many duplications as the second most supported location, which is in the root. Additionally, we observed an increase in duplications at the leaves of the species tree as p increased. Since most of the leaves are fixed episodes, the algorithm often assigned labels to create duplications at the leaves, resulting in larger sizes of episodes at leaves.
Conclusions and future outlook
In this article, we presented a novel problem that integrates genespecies mapping inference and genomic duplication detection. We proposed efficient algorithms to solve the problem exactly in the majority of instances, along with a heuristic modification for cases where exact solutions are not feasible. To demonstrate the effectiveness and accuracy of our proposed algorithm, we conducted computational experiments on both simulated and empirical data. While there is presently no established method for the simultaneous inference of duplications and genespecies mappings, the results showed that our algorithm was able to accurately infer recent WGD events when the number of missing labels was relatively small for simulated data. Moreover, the algorithm performed even better on empirical data, demonstrating its robustness and applicability to realworld scenarios. Nevertheless, our findings regarding genespecies mapping inference underscore the challenging nature of the problem with the current approach. Inferring true gene leaf labels proves difficult due to the combinatorial explosion of potential solutions and the resulting nearlyuniform distributions of genespecies mappings, which extend across substantial portions of the species tree.
To maximize topological similarities between a gene tree and its species tree, speciation nodes should more frequently appear in the resulting extensions of input partial gene trees. We observe that the optimization model tends to reconstruct leaf labels in a way that prioritizes duplication events assigned to the nearest fixed episodes or the root, in the absence of such episodes. This is confirmed by the property that fixed episodes are tight approximations of the EC cost, leading to a reduction in the number of speciation events in the final gene tree extensions. As a consequence, the model’s effectiveness may be limited in some cases when the number of unknown labels in partial gene trees is significant.
Several avenues for future research exist. For instance, the simulation pipeline could be extended to model fractionation [36, 48], offering a more detailed understanding of the biological processes involved. Alternatively, tools for segmental duplications [49] might provide a more accurate representation of gene duplication processes. Also, we plan to extend the analyzed model to strengthen the importance of the topological similarities between gene and species trees. Alternatively, one may limit the distance between the lcamapping of a gene duplication and its destination mapping in the final scenario similarly to [25]. Additionally, there are models of genomic duplications providing a higher level of detail than \(\textsf{EC}\), such as minimum episodes (ME) [20] and RMP [23], which can be adapted in a similar way to infer genespecies mappings and minimize the number of duplication episodes simultaneously. These models can be further combined with more general models of valid mappings, which allow the introduction of more duplication events than the minimum obtained by the lcamapping [4]. The combination of these models can provide a more comprehensive approach to inferring genespecies mappings and identifying the minimum number of duplication episodes.
Availability of data and materials
The software package metaEC, which is partially based on the embretnet repository, has been written in Python and is available with all datasets, including the commands to generate them and the output files, at https://bitbucket.org/pgor17/metaEC.
Notes
Here, the notation F[D] denotes the image of D on the function F.
Abbreviations
 EC:

Episode clustering
 DP:

Dynamic programming
 WGD:

Whole genome duplication
 CV:

Coefficient of variation
 MSA:

Multiple sequence alignment
 ME:

Minimum episodes
References
Goodman M, Czelusniak J, Moore GW, RomeroHerrera AE, Matsuda G. Fitting the gene lineage into its species lineage, a parsimony strategy illustrated by cladograms constructed from globin sequences. Syst Zool. 1979;28(2):132–63.
Page RDM. Maps between trees and cladistic analysis of historical associations among genes, organisms, and areas. Syst Biol. 1994;43(1):58–77.
Ma B, Li M, Zhang L. From gene trees to species trees. SIAM J Comput. 2000;30(3):729–52.
Górecki P, Tiuryn J. DLStrees: a model of evolutionary scenarios. Theoret Comput Sci. 2006;359(1–3):378–99.
Kuzmin E, VanderSluis B, Ba ANN, Wang W, Koch EN, Usaj M, Khmelinskii A, Usaj MM, Leeuwen J, Kraus O, Tresenrider A, Pryszlak M, Hu MC, Varriano B, Costanzo M, Knop M, Moses A, Myers CL, Andrews BJ, Boone C. Exploring wholegenome duplicate gene retention with complex genetic interaction analysis. Science. 2020;368(6498):5667.
Ohno S. Evolution by gene duplication. Berlin: Springer; 1970.
SalmanMinkov A, Sabath N, Mayrose I. Wholegenome duplication as a key factor in crop domestication. Nat Plants. 2016;2:16115.
Wu S, Lau KH, Cao Q, Hamilton JP, Sun H, Zhou C, Eserman L, Gemenet DC, Olukolu BA, Wang H, Crisovan E, Godden GT, Jiao C, Wang X, Kitavi M, ManriqueCarpintero N, Vaillancourt B, WiegertRininger K, Yang X, Bao K, Schaff J, Kreuze J, Gruneberg W, Khan A, Ghislain M, Ma D, Jiang J, Mwanga ROM, LeebensMack J, Coin LJM, Yencho GC, Buell CR, Fei Z. Genome sequences of two diploid wild relatives of cultivated sweetpotato reveal targets for genetic improvement. Nat Commun. 2018;9(1):4580.
Wolfe KH, Shields DC. Molecular evidence for an ancient duplication of the entire yeast genome. Nature. 1997;387(6634):708–13.
López S, Lim EL, Horswell S, Haase K, Huebner A, Dietzen M, Mourikis TP, Watkins TBK, Rowan A, Dewhurst SM, Birkbak NJ, Wilson GA, Loo PV, JamalHanjani M, Consortium T, Swanton C, McGranahan N. Interplay between wholegenome doubling and the accumulation of deleterious alterations in cancer evolution. Nat Genet. 2020;52(3):283–93.
Bielski CM, Zehir A, Penson AV, Donoghue MTA, Chatila W, Armenia J, Chang MT, Schram AM, Jonsson P, Bandlamudi C, Razavi P, Iyer G, Robson ME, Stadler ZK, Schultz N, Baselga J, Solit DB, Hyman DM, Berger MF, Taylor BS. Genome doubling shapes the evolution and prognosis of advanced cancers. Nat Genet. 2018;50(8):1189–95.
Quinton RJ, DiDomizio A, Vittoria MA, Kotýnková K, Ticas CJ, Patel S, Koga Y, Vakhshoorzadeh J, Hermance N, Kuroda TS, Parulekar N, Taylor AM, Manning AL, Campbell JD, Ganem NJ. Wholegenome doubling confers unique genetic vulnerabilities on tumour cells. Nature. 2021;590(7846):492–7.
Guigó R, Muchnik IB, Smith TF. Reconstruction of ancient molecular phylogeny. Mol Phylogenet Evol. 1996;6(2):189–213.
Page RDM, Cotton JA. Vertebrate phylogenomics: reconciled trees and gene duplications. Pacific Symposium on Biocomputing. 2002:536–47.
Fellows M, Hallet M, Stege U. On the multiple gene duplication problem. In: 9th International Symposium on Algorithms and Computation (ISAAC’98), Lecture Notes in Computer Science 1533, Taejon, Korea, 1998:347–356.
Bansal MS, Eulenstein O. The multiple gene duplication problem revisited. Bioinformatics. 2008;24(13):132–8.
Burleigh JG, Bansal MS, Wehe A, Eulenstein O. Locating multiple gene duplications through reconciled trees. In: Research in Computational Molecular Biology: 12th Annual International Conference, RECOMB 2008, Singapore, March 30April 2, 2008. Proceedings 12, 2008:273–284. Springer
Luo CW, Chen MC, Chen YC, Yang RWL, Liu HF, Chao KM. Lineartime algorithms for the multiple gene duplication problems. IEEE/ACM Trans Comput Biol Bioinf. 2011;8(1):260–5.
Mettanant V, Fakcharoenphol J. A lineartime algorithm for the multiple gene duplication problem. In: The 12th National Computer Science and Engineering Conference (NCSEC), 2008:198–203.
Paszek J, Górecki P. Efficient algorithms for genomic duplication models. IEEE/ACM Trans Comput Biol Bioinf. 2018;15(5):1515–24.
Paszek J, Górecki P. Genomic duplication problems for unrooted gene trees. BMC Genomics. 2016;17(1):165–75.
Paszek J, Górecki P. Inferring duplication episodes from unrooted gene trees. BMC Genomics. 2018;19(S5).
Iersel LV, Janssen R, Jones M, Murakami Y, Zeh N. PolynomialTime Algorithms for Phylogenetic Inference Problems involving duplication and reticulation. IEEE/ACM Trans Comput Biol Bioinf. 2019
Paszek J, Tiuryn J, Górecki P. Minimizing genomic duplication episodes. Comput Biol Chem. 2020;89: 107260.
Paszek J, Markin A, Górecki P, Eulenstein O. Taming the duplicationlosscoalescence model with integer linear programming. J Comput Biol. 2021;28(8):758–73.
Dondi R, Lafond M, Scornavacca C. Reconciling multiple genes trees via segmental duplications and losses. Algorithms Mol Biol. 2019;14:7.
RoyoLlonch M, Sánchez P, RuizGonzález C, Salazar G, PedrósAlió C, Sebastián M, Labadie K, Paoli L, Ibarbalz FM, Zinger L, Churcheward B, Coordinators TO, Chaffron S, Eveillard D, Karsenti E, Sunagawa S, Wincker P, KarpBoss L, Bowler C, Acinas SG. Compendium of 530 metagenomeassembled bacterial and archaeal genomes from the polar Arctic Ocean. Nat Microbiol. 2021;6(12):1561–74.
Wirbel J, Pyl PT, Kartal E, Zych K, Kashani A, Milanese A, Fleck JS, Voigt AY, Palleja A, Ponnudurai R, Sunagawa S, Coelho LP, SchrotzKing P, Vogtmann E, Habermann N, Niméus E, Thomas AM, Manghi P, Gandini S, Serrano D, Mizutani S, Shiroma H, Shiba S, Shibata T, Yachida S, Yamada T, Waldron L, Naccarati A, Segata N, Sinha R, Ulrich CM, Brenner H, Arumugam M, Bork P, Zeller G. Metaanalysis of fecal metagenomes reveals global microbial signatures that are specific for colorectal cancer. Nat Med. 2019;25(4):679–89.
Betkier A, Szczęsny P, Górecki P. Fast algorithms for inferring genespecies associations. In: Bioinformatics Research and Applications: 11th International Symposium, ISBRA 2015 Norfolk, USA, June 7–10, 2015 Proceedings 11, 2015:36–47. Springer.
Zhang L, Cui Y. An efficient method for dnabased species assignment via gene tree and species tree reconciliation. In: Algorithms in Bioinformatics: 10th International Workshop, WABI 2010, Liverpool, UK, September 6–8, 2010. Proceedings 10, 2010:300–311. Springer.
Mykowiecka A, Szczęsny P, Górecki P. Inferring genespecies assignments in the presence of horizontal gene transfer. IEEE/ACM Trans Comput Biol Bioinf. 2017;15(5):1571–8.
Łukasiewicz J. Selected Works, vol. 1. Amsterdam: NorthHolland Publishing Company; 1970.
Mallo D, De Oliveira Martins L, Posada D. Simphy: phylogenomic simulation of gene, locus, and species trees. Syst Biol. 2016;65(2):334–44.
Molloy EK, Warnow T. FastMulRFS: fast and accurate species tree estimation under generic gene duplication and loss models. Bioinformatics. 2020;36(Supplement1):57–65.
Rasmussen MD, Kellis M. Unified modeling of gene duplication, loss, and coalescence using a locus tree. Genome Res. 2012;22(4):755–65.
Cheng F, Wu J, Cai X, Liang J, Freeling M, Wang X. Gene retention, fractionation and subgenome differences in polyploid plants. Nat Plants. 2018;4:258–68.
Fletcher W, Yang Z. Indelible: a flexible simulator of biological sequence evolution. Mol Biol Evol. 2009;26(8):1879–88.
Price MN, Dehal PS, Arkin AP. FastTree: computing large minimum evolution trees with profiles instead of a distance matrix. Mol Biol Evol. 2009;26(7):1641–50.
Górecki P, Tiuryn J. Urec: a system for unrooted reconciliation. Bioinformatics. 2007;23(4):511–2.
Górecki P, Rutecka N, Mykowiecka A, Paszek J. Simultaneous Reconstruction of Duplication Episodes and GeneSpecies Mappings. In: Belazzougui D, Ouangraoua A, editors. 23rd International Workshop on Algorithms in Bioinformatics (WABI 2023), vol. 273. Leibniz International Proceedings in Informatics (LIPIcs). Dagstuhl, Germany: Schloss Dagstuhl—LeibnizZentrum für Informatik; 2023. p. 6–1618.
Yilmaz B, Elibol E, Shangpliang HNJ, Ozogul F, Tamang JP. Microbial communities in homemade and commercial kefir and their hypoglycemic properties. Fermentation. 2022;8(11):590.
Feng B, Lin Y, Zhou L, Guo Y, Friedman R, Xia R, Hu F, Liu C, Tang J. Reconstructing yeasts phylogenies and ancestors from whole genome data. Sci Rep. 2017;7(1):1–12.
MarcetHouben M, Gabaldón T. Beyond the wholegenome duplication: phylogenetic evidence for an ancient interspecies hybridization in the baker’s yeast lineage. PLoS Biol. 2015;13(8):1002220.
Consortium TU. Uniprot: the universal protein knowledgebase in 2023. Nucleic Acids Res. 2023;51(D1):523–31.
Van Dongen S. Graph clustering via a discrete uncoupling process. SIAM J Matrix Anal Appl. 2008;30(1):121–41.
Edgar RC. Muscle: a multiple sequence alignment method with reduced time and space complexity. BMC Bioinf. 2004;5(1):1–19.
Guindon S, Dufayard JF, Vincent L, Anisimova M, Hordijk W, Gascuel O. New algorithms and methods to estimate maximumlikelihood phylogenies: assessing the performance of phyml 3.0. Syst Biol. 2010;59(3):307–21.
Zhang Y, Zheng C, Sankoff D. Pinning down ploidy in paleopolyploid plants. BMC Genomics. 2018;19(Suppl 5):287.
Davín AA, Tricou T, Tannier E, Vienne DM, Szöllősi GJ. Zombi: a phylogenetic simulator of trees, genomes and sequences that accounts for dead linages. Bioinformatics. 2020;36(4):1286–8.
Acknowledgements
The authors would like to thank the anonymous reviewers for their valuable suggestions and comments.
Funding
The support was provided by National Science Centre grant nr #2017/27/B/ST6/02720. JP was supported by National Science Centre Grant #2020/39/D/ST6/03321. Publication costs were covered by the funds of National Science Centre Grant #2020/39/D/ST6/03321.
Author information
Authors and Affiliations
Contributions
PG developed and implemented the algorithms. NR and PG performed the simulations and computational tests. AM performed the empirical evaluation. NR contributed to the development of the software. JP contributed to the simulated dataset WGD analysis and development of algorithms, and all authors contributed to the writing. All authors reviewed and approved the manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Górecki, P., Rutecka, N., Mykowiecka, A. et al. Unifying duplication episode clustering and genespecies mapping inference. Algorithms Mol Biol 19, 7 (2024). https://doi.org/10.1186/s13015024002528
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s13015024002528