Algorithms of Ancestral Gene Length Reconstruction

Abstract

Ancestral sequence reconstruction is a well-known problem in molecular evolution. The problem presented in this study is inspired by sequence reconstruction, but instead of leaf-associated sequences we consider only their lengths. We call this problem ancestral gene length reconstruction. It is a problem of finding an optimal labeling which minimizes the total length’s sum of the edges, where both a tree and nonnegative integers associated with corresponding leaves of the tree are the input. In this paper we give a linear algorithm to solve the problem on binary trees for the Manhattan cost function 𝑠(𝑣,𝑤)=|𝜋(𝑣)−𝜋(𝑤)|.

1. Introduction

Ancestral sequence reconstruction (ASR) is a well-recognized problem in molecular evolution [1]. Let 𝐆 be a (phylogenetic) tree with 𝐧 leaf nodes, and 𝑘 strings over one alphabet (gene sequences) assigned to 𝑘 leaves (𝑘≤𝑛). ASR may be defined in the following way: assignment of strings to inner nodes “in the best possible way.” There are two main paradigms for ASR: maximum parsimony (MP) and probabilistic-based reconstruction. The latter includes maximum likelihood (ML) and Bayesian reconstructions. MP reconstruction has a time complexity linear in the number of sequences analyzed. The problem of the parsimonious reconstruction of ancestral states for the given tree with the given states of its leaves (the most parsimonious assignment of the labels of internal nodes for a fixed tree topology) is a well-studied problem [2–4]. Efficient algorithms have also been developed for different types of ML-based reconstructions (reviewed in [5]). ASR methods require as input both a phylogenetic tree and a set of gene sequences associated with corresponding leaves of the tree [6].

ASR is related to gene sequence evolution while the problem presented in this paper, being inspired by ASR, deals with gene length variation. Instead of considering leaf-associated sequences we take into account only their lengths. Instead of the reconstruction of ancestral sequences, we search for the optimal reconstruction of ancestral gene lengths. The problem may be called ancestral gene length reconstruction (AGLR). AGLR is actually a problem of finding an optimal labeling which minimizes the total “length” sum of the edges, the minimum sum problem where both a tree and nonnegative integers associated with corresponding leaves of the tree are the input.

In the graph theory vertex labeling related problems were intensively studied [8]. Typically, the problems can be described as follows: for a given graph, find the optimal way of labeling the vertices with distinct integers. The problems and their solutions were described in [9–12]. In [13] we presented the algorithms to solve the minimum sum problem where both a tree and positive integers associated with all leaves of the tree are the input (finding the optimal way of labeling the vertices with positive integers). Here we would like to formulate the minimum sum problem where both a tree and positive integers associated with some of the leaves of the tree are the input (finding the optimal way of labeling the vertices with nonnegative integers). This problem reflects a situation in which the genome tree is constructed by one or another method for a set of genomes, the leaves of the tree are linked with the corresponding genomes of the set, and the leaves are labeled by integers designating lengths of genes of a chosen gene family. Some leaves would be labeled zero because corresponding genomes have no genes of the chosen gene family. Alternatively, it may be a case of a missing value but in this study we do not consider this case: in the problem definition that we bring here zero means “no value.”

In this paper we provide a linear algorithm to solve max sum problem on binary trees for the Manhattan cost function 𝑠(𝑣,𝑤)=|𝜋(𝑣)−𝜋(𝑤)|. The algorithm uses dynamic programming technique and the properties of the Manhattan distance.

2. Preliminaries

Let 𝐆 be a tree with 𝐧 leaf nodes, vertex set 𝐕(𝐆), and edge set 𝐄(𝐆). 𝑁=|𝑉(𝐺)|. Let us number the leaf nodes of 𝐺:1,2,…,𝑛. Let us number the root of 𝐺∶𝑁. An integer labeling 𝝅 of 𝐺 is a mapping 𝜋 from 𝐺 to a set of nonnegative integers, where label 𝟎 is an out-of-the ordinary label meaning “absent value.” Let us denote integer labeling of the leaf nodes of 𝐺(𝜋(1)=𝑝₁,…,𝜋(𝑛)=𝑝_𝑛). Let us denote by 𝑔_min and 𝑔_max minimum and maximum positive integers labeling leaf nodes: 𝑔_min=min𝑝_𝑖∶𝑝_𝑖>0; 𝑔_max=max𝑝_𝑖; 𝑚=𝑔_max−𝑔_min+1.

Let us introduce a cost function 𝜑 of the edge 𝑣𝑤∈𝐸(𝐺): 𝜑=𝐶𝑥,𝑦0if𝑥=𝑦else₁𝐶if𝑥=0else₂𝜃,if𝑦=0else𝑥,𝑦 (1) where the nonnegative cost function 𝜃(𝑥,𝑦) has the following distance properties:(i)𝜃𝑥,𝑦≥0𝑥=𝑦⟷𝜃𝑥,𝑦=0/∗functionisequaltozeroifandonlyifitsargumentsareequal∗/ (2)(ii)𝜃𝑥,𝑦=𝜃𝑦,𝑥/∗symmetry∗/ (3)(iii)𝜃,𝜃;𝑥>𝑦⟶𝑥,𝑦<𝜃𝑥,𝑦−1𝑥,𝑦<𝜃𝑥+1,𝑦 (4)(iv)𝜃,𝜃.𝑥<𝑦⟶𝑥,𝑦<𝜃𝑥−1,𝑦𝑥,𝑦<𝜃𝑥,𝑦+1 (5)𝐶₁>𝐶₂>𝑚=𝑔_max−𝑔_min+1. 𝐶₁ is a gain penalty, 𝐶₂ is a loss penalty, and 𝜃 is a length change penalty function. Since the likelihoods of loss and gain events are likely to differ, we may need to weight them differently. This is achieved by introducing different penalties 𝐶₁>𝐶₂; the loss penalty is normally assigned a value close to 𝑔_max−𝑔_min, whereas the gain penalty should be larger due to biological considerations. They suggest that, on average, gene loss might be a more likely event than gene gain. Therefore, different gain penalties were used in our study similarly to as it was done in [14].

An example of a function 𝜃(𝑥,𝑦) is |𝜋(𝑣)−𝜋(𝑤)|^𝜆. In case of 𝜆=1 we obtain an absolute value of the difference between labelings 𝑣 and 𝑤: |𝜋(𝑣)−𝜋(𝑤)|. In case of 𝜆=2 we obtain a square of the difference between labelings 𝑣 and 𝑤: (𝜋(𝑣)−𝜋(𝑤))².

2.1. An Arbitrary Tree and an Arbitrary Cost Function

Given a tree 𝐺, an integer labeling of the leaves of 𝐺(𝑝₁,…,𝑝_𝑛)=1, the gain penalty 𝐶₁, the loss penalty 𝐶₂, and a cost function 𝜃 ((1)–(5)), the minimum sum problem is to find a labeling which minimizes the total cost: 𝑆(𝐺)=min_{𝜋(𝐺)∀{𝑣𝑤}∈𝐸}𝜑(𝜋(𝑣),𝜋(𝑤))overall𝜋. (6)

2.2. A Binary Tree Problem

Given a binary tree 𝐺, an integer labeling of the leaves of 𝐺(𝑝₁,…,𝑝_𝑛), the “gain” penalty 𝐶₁, and the “loss” penalty 𝐶₂, the Manhattan minimum sum problem is to find the labelings which minimize the sum 𝑆 over all 𝜋𝑆(𝐺)=_{∀{𝑣𝑤}∈𝐸(𝐺)&𝜋(𝑣)≠0&𝜋(𝑤)≠0}|𝜋(𝑣)−𝜋(𝑤)|+𝑘₁⋅𝐶₁+𝑘₂⋅𝐶₂, (7) where 𝑘₁ is a number of edges of type (𝜋(𝑣)=0&𝜋(𝑤)>0), and 𝑘₂ is a number of edges of type (𝜋(𝑣)>0&𝜋(𝑤)=0).

3. Problem Solutions

3.1. DP Algorithm (for the Problem (1))

Due to the properties ((2)–(5)) of the cost function 𝜃(𝑥,𝑦) all labels of the optimal labeling must be either equal to 0 or in the interval [𝑔_min,𝑔_max]. As a consequence of this, the dynamic programming (DP) method is applicable for the problem. It will be easier to explain the DP method on a binary tree using 𝜎_𝑘(𝑖) notation. The quantity 𝜎_𝑘(𝑖) will be interpreted as the minimal cost, given that node 𝑘 is assigned integer 𝑖, to the subtree with the node 𝑘 as a root of the subtree.

3.1.1. DP Algorithm for a Binary Tree

Up Phase. A procedure called DP_up calculates the costs 𝜎_𝑘(𝑖) of all nodes 𝑉(𝐺) of the tree 𝐺, given a cost function φ.

When we compute 𝜎_𝑁(𝑖) for the root node (the index of the root is 𝑁), then we simply choose the minimum of these values: 𝑆(𝐺)=min_𝑖𝜎_𝑁(𝑖) (8)

Initiation. Given labeling of the leaf nodes of 𝐺(𝑝₁,…,𝑝_𝑛} at the tips of the tree the 𝜎_𝑖(𝑗) are easy to compute. The cost is 0 if the observed integer 𝑝_𝑖 is integer 𝑗, and infinite otherwise. 𝜎_𝑖𝑗=0if𝑗=𝑝_𝑖∞otherwise. (9)

Iteration. For the immediate common ancestor of the nodes 𝑙 and 𝑟, node 𝑎, we have 𝜎_𝑎(0)=min𝜎_𝑙(0),𝐶₁+min_𝑗𝜎_𝑙𝑗+min𝜎_𝑟(0),𝐶₁+min_𝑘𝜎_𝑟(𝑘),𝜎_𝑎(𝑖)=minmin_𝑗𝜃𝑖,𝑗+𝜎_𝑙𝑗,𝐶₂+𝜎_𝑙(0)+minmin_𝑘𝜃(𝑖,𝑘)+𝜎_𝑟(𝑘),𝐶₂+𝜎_𝑟(0),∀𝑖,𝑗,𝑘∈𝑔_min,𝑔_max. (10)

The interpretation of this equation is immediate. The smallest possible cost given that node 𝑎 is assigned zero is either the cost 𝜎_𝑙(0) or the “gain” penalty 𝐶₁ plus the minimum of 𝜎_𝑙(𝑗), the least of the two plus the minima of corresponding values associated with the right descendant tree. The smallest possible cost given that node 𝑎 is assigned 𝑖 is a sum of two values: the first one is either the cost 𝜃(𝑥,𝑦) of the edge from node 𝑎 to node 𝑙, plus the cost 𝜎_𝑙(𝑗) of the left descendant subtree given that node 𝑙 is in state 𝑗, or the “loss” penalty 𝐶₂ plus 𝑆_𝑙(0); the second one is the cost 𝜃(𝑖,𝑘) of the edge from the node 𝑎 to the node 𝑟, plus the cost 𝜎_𝑟(𝑘) of the right descendant subtree given that node 𝑟 is in state 𝑘. We select those values of 𝑗 and 𝑘 which minimize that sum. Equation (10) is applied successively to each inner node in the tree, doing a postorder tree traversal. Finally it computes all the 𝜎_𝑁(𝑖), and then (8) is used to find the minimum cost for the whole tree. The complexity of the Up_phase of the algorithm is 𝑂(𝑁^∗𝑚^∗𝑚).

Traceback. The procedure calculates the labels 𝜋(𝑝) of all nodes 𝑝 of the tree 𝐺.

Choose any integer 𝑖 which provides the minimum of the 𝜎_𝑁(𝑖)—it is the root label. It may be either zero or a positive 𝑖. Doing a preorder tree traversal, successively label each inner node in the tree: for any inner node 𝑝, and given that a parent label 𝑖 was reconstructed, the label 𝜋(𝑝)=𝑗 is easily reconstructed as well.

3.1.2. DP Algorithm for an Arbitrary Tree

Up-Phase. A procedure DP_up calculates the costs 𝜎_𝑘(𝑖) of all nodes 𝑉(𝐺) of the tree.

Suppose that the 𝑘_𝑎 descendant nodes of the node 𝑎 are called 𝑏_𝑗. The following equation will therefore be similar to (10) replacing the sum of 𝜎_𝑙 and 𝜎_𝑟 by the total sum of 𝜎_𝑗1, while 𝑗₁ traverses all values of 𝑏_𝑗: 𝜎_𝑎(0)=^𝑘𝑎_𝑗1min𝜎_𝑗1(0),𝐶₁+min_𝑗𝜎_𝑗1𝑗,𝜎_𝑎(𝑖)=^𝑘𝑎_𝑗1minmin_𝑗𝜃𝑖,𝑗+𝜎_𝑗1𝑗,𝐶₂+𝜎_𝑗1(0,𝑔)∀𝑖,𝑗∈_min,𝑔_max. (11)(12)

This equation is applied successively to each node in the tree, doing a postorder tree traversal. Finally it computes all the 𝜎_𝑁(𝑖), and then (8) is used to find the minimum cost for the whole tree.

Down Phase. As Traceback above: Consider the following.

3.2. DP Algorithm for a Manhattan Sum for a Binary Tree (Problem (2))

Manhattan distance 𝜃(𝜋(𝑣),𝜋(𝑤)) is an absolute value of the difference between labelings 𝑣 and 𝑤∶|𝜋(𝑣)−𝜋(𝑤)|. This distance measure has the following property: if siblings have positive labels, then all integers that lie between these values may equally serve as optimal labels of a parent.(i)If (𝜋(𝑙)≤𝜋(𝑟)), then for all 𝑘𝜋(𝑙)≤𝑘≤𝜋(𝑟) the score 𝜃(𝑘,𝜋(𝑙))+𝜃(𝑘,𝜋(𝑟))=𝑘−𝜋(𝑙)+𝜋(𝑟)−𝑘=𝜋(𝑟)−𝜋(𝑙).(ii)If (𝜋(𝑙)≤𝜋(𝑟)), then for all 𝑘<𝜋(𝑙)≤𝜋(𝑟) the score 𝜃(𝑘,𝜋(𝑙))+𝜃(𝑘,𝜋(𝑟))=𝜋(𝑙)−𝑘+𝜋(𝑟)−𝑘=𝜋(𝑟)−𝜋(𝑙)+2(𝜋(𝑙)−𝑘).(iii)If (𝜋(𝑙)≤𝜋(𝑟)), then for all 𝜋(𝑙)≤𝜋(𝑟)<𝑘 the score 𝜃(𝑘,𝜋(𝑙))+𝜃(𝑘,𝜋(𝑟))=𝑘−𝜋(𝑙)+𝑘−𝜋(𝑟)=𝜋(𝑟)−𝜋(𝑙)+2(𝑘−𝜋(𝑟)).

So, as it would be proven below, at the bottom-up stage of the DP algorithm it would be sufficient to assign to each node 𝑎 in the tree 𝐺 four values: left(𝑎), right(𝑎), 𝑍(𝑎), and 𝑋(𝑎). The meanings of the values are as follows: left and right are bounds of an interval associated with the node 𝑎, 𝑍 is a cost value 𝜎_𝑎(0), and 𝑋 is a cost 𝜎_𝑎(𝑖) for any integer 𝑖 from the interval: left ≤ 𝑖 ≤ right.

Initiation. Given labeling of the leaf nodes of 𝐺(𝑝₁,…,𝑝_𝑛)=1 these four values are easy to compute for the leaf nodes: for (𝑖=1;𝑖≤𝑛;𝑖++) if (𝑝[𝑖]==0) {𝑍[𝑖]=0; left[𝑖]=0; right[𝑖]=0; 𝑋[𝑖]=𝐶₁+𝐶₂} else {𝑍[𝑖]=𝐶₁+𝐶₂; left[𝑖]=𝑝[𝑖]; right[𝑖]=𝑝[𝑖];𝑋[𝑖]=0}.

3.2.1. Examples

Let us consider the simplest trees with two, three, and four labeled leaves. The simplest tree configuration is presented in Figure 1. There is only one node to label—the root node.(i)Figure 1(a): no genes are assigned to the leaves → no gene is assigned to the root.(ii)Figure 1(b): the left leaf has no gene, and the right leaf has a gene with the length equal to 136 → the root is labeled by 136; the score is equal to the loss penalty 𝐶₂=30.(iii)Figure 1(c): any label 125≤𝑘≤136 is good to label the root; the score is equal to 136−125=11.

(a)

(b)

(c)

Figure 1 

Assignment of bottom-up stage values (left, right, 𝑍, and 𝑋) in 2-leaf trees. The “gain” penalty 𝐶1=50; the “loss” penalty 𝐶2=30. Optimal labels are in red.

The next simplest tree topology—three-leaf trees—is presented in Figure 2. There are two nodes to label, the inner node and the root.(i)Figure 2(a): the inner node is labeled analogically to the root in Figure 1(c): any 𝑘  125≤𝑘≤136 is equally good to label the inner node; the root node is labeled analogically to the root in Figure 1(b): (𝑍(root)=𝐶₁+(136−125))>(𝑋(root)=𝐶₂+11) → the root is labeled by any 𝑘  125≤𝑘≤136, that is, by 125.(ii)Figure 2(b): labeling is similar to that of Figure 1(a).(iii)Figure 2(c): the inner node is labeled analogically to Figure 2(a): any label 125≤𝑘≤136 is good to label it; the score is equal to 136−125=11. The root should be labeled by 136 because 125<136<141.

(a)

(b)

(c)

Figure 2 

Assignment of bottom-up stage values (left, right, 𝑍, and 𝑋) in 3-leaf trees. The “gain” penalty 𝐶1=50; the “loss” penalty 𝐶2=30. Optimal labels are in red.

Determination of the optimal labeling of the four-leaf trees is very similar to the examples described above. Figure 3 illustrates labeling of the tree where all four leaves have nonzero labels: ((125, 141), (136, 150)). Labeling of the inner nodes is as above (Figure 2(c)): [125,141] and [136,150]. All integers of the intersection between these two close intervals are optimal values to label the root: [125,141]∩[136,150]=[136,141]. In Figure 3 we present the value 136 as a chosen suitable label.

Figure 3 

Assignment of bottom-up stage values (left, right, 𝑍, and 𝑋) in a 4-leaf tree with all four leaves labeled by positive integers. The “gain” penalty 𝐶1=50; the “loss” penalty 𝐶2=30. Optimal labels are in red.

Examples of the trees with very distinct subtrees are presented in Figures 4 and 5. In Figure 4 we present a tree obtained by merging two very different subtrees. The left 4-leaf subtree has very obvious intuitive labeling of internal nodes: all nodes should be labeled by zero. The right subtree is identical to the tree presented in Figure 2(c). Merging of these two subtrees produces bottom-up stage values (left, right, 𝑍, and 𝑋) to the new root equal to [125,136,111,91]. In spite of assignins the interval [125,136] to the root only the value 136 provides the optimal solution. (We would like to express our gratitude to the anonymous reviewer for bringing our attention to this situation.) We formulate this rule below describing traceback stage of the algorithm. Figure 4 is chosen to illustrate labeling of nodes similar to the root of the tree.

Figure 4 

Labeling of a “peculiar” tree. The left subtree has three zero and one nonzero leaf, while the right subtree has three nonzero leaves. The “gain” penalty 𝐶1=50; the "loss" penalty 𝐶2=30. Optimal labels are in red.

Figure 5 

Archaeal part of Figure 4(b) from [7] labeled accordingly to COG0835.

After considering these few simple examples, we describe the algorithm.

3.2.2. Bottom-Up Stage

Initiation. Given labeling of the leaf nodes of 𝐺(𝑝₁,…,𝑝_𝑛) at the tips of the tree the 𝜎_𝑖(𝑗) are easy to compute. The cost is 0 if the observed integer 𝑝_𝑖 is integer 𝑗, and 𝐶₁+𝐶₂otherwise. (13)

Iteration. Doing a postorder tree traversal assign successively to each node in the tree the abovementioned four values left(𝑎), right(𝑎), 𝑍(𝑎), and 𝑋(𝑎). An interval [left(𝑎), right(𝑎)] is assigned according to the following rule: if anyone of two children intervals is not defined, then assign the interval of the other child; otherwise, a parent interval is either an intersection of the intervals of its children or an interval that lies between these intervals if their intersection is empty. 𝑍 is a cost value 𝜎_𝑎(0), where for the Manhattan distance we can rewrite (10) as 𝑍(𝑎)=𝜎_𝑎(0)=min𝜎_𝑙(0),𝐶₁+𝜎_𝑙𝑗+min𝜎_𝑟(0),𝐶₁+𝜎_𝑟𝑗,𝑋(𝑎)=𝜎_𝑎(𝑖)=minmin_𝑗𝜃𝑖,𝑗+𝜎_𝑙𝑗,𝐶₂+𝜎_𝑙(0)+minmin_𝑘𝜃(𝑖,𝑘)+𝜎_𝑟(𝑘),𝐶₂+𝜎_𝑟(0)=minmin_𝑗𝑖−𝑗+𝜎_𝑙𝑗,𝐶₂+𝜎_𝑙(0)+minmin_𝑗𝑖−𝑗+𝜎_𝑟𝑗,𝐶₂+𝜎_𝑟(0). (14)

3.2.3. Pseudocode

For more details see Pseudocode 1.

/* Assignment values left[𝑎] and right [𝑎] */ FL = FALSE; if (left[l] == 0) { left[a] = left[r]; right[a] = right[r]} else if (left[r] == 0) { left[a] = left[l]; right[a] = right[l]} else if (left[l] ≤ left[r]) {   if (right[l] < left[r]) {    left[a] = right[l]; right[a] = left[r]; FL = TRUE;   } else {    left[a] = left[r]; right[a] = min(right[l], right[r]);   } } else {  if (right[r] < left[l]) {    left[a] = right[r]; right[a] = left[l]; FL = TRUE;   } else {    left[a] = left[l]; right[a] = min(right[l], right[r]);   } } /* Assignment values 𝑍[𝑎] and 𝑋[𝑎]: 𝑍 is a cost of 𝜎𝑎(0), 𝑋 is a cost of 𝜎𝑎(𝑖) for all i: left ≤ 𝑖 ≤ right */ if (FL) diff = right[𝑎] − left[𝑎]; else diff = 0; 𝑍(𝑎) = min(𝑍(𝑙), 𝐶1+𝑋(𝑙)) + min(𝑍(𝑟),𝐶1+𝑋(𝑟)) 𝑋(𝑎) = min(diff + 𝑋(𝑙) + 𝑋(𝑟), 𝑋(𝑟) + 𝐶2+ 𝑍(𝑙)), 𝑋(𝑙) + 𝐶2+ 𝑍(𝑟), 2𝐶2 + 𝑍(𝑙) + 𝑍(𝑟))

Pseudocode 1 

3.2.4. Traceback Stage

Interval Correction Rule. Following the bottom-up stage four values left(𝑎), right(𝑎), 𝑍(𝑎), and 𝑋(𝑎) are assigned to every internal node 𝑎 of the tree. An interval (left(𝑎), right(𝑎)) should be diminished if one of the edges connecting the node 𝑎 with its son becomes of type (𝑘,0), 𝑘>0. Let us denote sons of the node 𝑎 by 𝑙(𝑎) and 𝑟(𝑎). Correction condition Ω(𝑎) would be formulated as Ω(𝑎)=(𝑋(𝑎)(𝑎𝑋(𝑙(𝑎(𝑙(𝑎𝑉(𝑋(𝑟(𝑎(𝑟(𝑎.≤𝑍))and(())>𝑍)))))>𝑍)))) (15) If Ω(𝑎) is TRUE, then the bounds of the corrected interval would be obtained by intersection of the interval associated with the node with the corrected interval associated with the corresponding son: if𝑋(𝑙(𝑎))>𝑍(𝑙(𝑎)),thenleft^′(𝑎),right^′(𝑎=(𝑎)(𝑎∩(𝑟(𝑎(𝑟(𝑎)left,right)left)),right))elseleft^′(𝑎),right^′(𝑎=(𝑎)(𝑎∩(𝑙(𝑎(𝑙(𝑎;)left,right)left)),right)) (16) Otherwise, the bounds of the corrected interval would not be changed from the original ones: left^′(𝑎),right^′(𝑎)=left(𝑎),right(𝑎). (17)

Initiation. Labeling of the root: if 𝑋(𝑁)≤𝑍(𝑁), then correct the root interval according to (15)–(17), and then choose an integer from the corrected interval assigned to the root node otherwise choose 0—it is the root label 𝜋(𝑁).

Iteration. Doing a preorder tree traversal, successively label each node in the tree either by an integer from the corrected interval assigned to this node which is the nearest to its parent label (it may be either the value equal to the parent label or the boundary value of the interval assigned with the node) or by 0.

The proof of the correctness of a simpler algorithm (without zero-labeled leaves) is published in [15]. In the Appendix there are several lemmas, from which the correctness of the algorithm presented here results.

3.2.5. Example

In Figure 5 of [7] the consensus trees obtained from 100 genome trees were presented. The trees were produced on the basis of 80% randomly chosen COGs, and the right tree was produced on the basis of 15%-jackknifing (the explanations in the text of [7]). This tree possesses phylogenetic reasonableness.(a)The representatives of both prokaryotic Kingdoms: Eubacteria and Archaea are clustered separately. In other words, Archaeal organisms (genomes 0, 1, 8, 29–32, and 35–50) form a monophyletic group.(b)Euryarchaeota and Crenarchaeota form monophyletic groups.

A part of this tree was selected to illustrate the algorithm. We took the upper part of the tree related exclusively to Archaea (see A/B marked arrow in Figure 4(b) from [7]) and placed the root at the point dividing all Archaeal genomes into Euryarchaeota and Crenarchaeota (see E/C marked arrow in Figure 4(b) from [7]). Thus, Figure 5 is a part of Figure 4(b) from [7] labeled according to COG0835. This COG was randomly selected as suitable for purposes of illustration. Table 1 presents a list of Archaeal genomes from the whole set of genomes that were used for a genome tree construction (Figure 4(b) from [7]). Table 2 presents the lengths of the Archaeal proteins of this COG.

Table 1 

List of archaeal genomes for Figure 4.

Table 2 

Protein lengths of the chemotaxis signal transduction proteins. Archaeal part of COG0835.

To assign labels to the leaves of the tree of Figure 5 two preprocessing steps were done: (1) taking off outliers, the lengths 328 of the 𝐻. marismortui protein and 344 of the M. hungatei protein are obvious outliers; (2) taking the median value of paralog’s lengths of the genomes 30, 31, 46, 48, and 50. Figure 5 presents results of application of the bottom-up and traceback stages of the algorithm to this tree: a quartet that was assigned to a node 𝑎 at the bottom upstage is shown under the edge linking the node 𝑎 and its parent node, a label that was assigned to the node 𝑎 at the traceback stage, is shown over the same edge.

As we can see the root is labeled by zero. There are two gene-birth events and one gene-death event. One gene was born with the length of 155 and another gene birth is labeled by 146. Genome number 32 (Haloquadratum walsbyi) has no protein from COG0835, while other Haloarchaea (genomes 29–31) do have. Thus, the edge connecting with leaf labeled by 32 is marked with a gene-loss symbol.

4. Discussion

In [15] the algorithms to find the optimal labeling of the vertices of the tree under Wagner parsimony were presented. A simple extension of the problem could be finding the optimal labeling of the vertices of the tree with nonnegative integers. This more realistic approach requests special consideration of zero labeling. Wedges of type (𝑘,0), 𝑘>0, should be scored differently from wedges of type (0,𝑘), 𝑘>0, because the (𝑘,0) notes gene loss, while (0,𝑘) notes gene gain. These events should be scored differently. Interestingly, this differentiated scoring in addition to tree labeling resulted in reconstruction of “parsimonious” evolutionary scenario. Reconstruction of a gene evolution along a species tree is an interesting and principal problem. Lyubetsky and his coworkers contributed a lot to formulation and solving this problem. In their studies [16–22] the authors tackled mainly two important and sophisticated phylogenetic problems. The obtained results are partially reviewed in the first section of [22] which also provides an extended biological background and relevant references. Reconstruction of a gene evolution along a species tree (to build the evolutionary scenario), following the approach of Lyubetsky et al., is to find an optimal mapping of a gene tree into a species tree. (An example of a different approach was presented in [14].) The second problem is to construct a supertree from the given set of gene trees.

As it was mentioned in [22], the first problem, stated as a tree-into-tree mapping, is solved in polynomial (often linear, and at maximum cubic) time even for the case of time slices and horizontal gene transfers. The algorithms presented in our study are polynomial as well.

Choosing 𝐶₁ (a gain penalty), 𝐶₂ (a loss penalty), and 𝜃 (a label change penalty) is crucial for reconstruction of trustworthy evolutionary scenario. However, it is very difficult task and we cannot claim categorically that choosing “correct” parameters of the model will result in truly reliable reconstruction. We do plan to make a comparison between results obtained by abovementioned methods of Lyubetsky and ours (work in progress).

To prepare input for the algorithm, as it was done above for 3.2.5, the original data is to be transformed to the following format: to each (genome, COG) pair one standardized protein length should be assigned (as we described in [7]). For a given COG, each organism is represented by a calculated length—a median length of all paralogous proteins. A natural extension would be to formulate the labeling problem taking into account existence of paralogs.

We may define a 𝑘-tuple integer labeling Π of 𝐺 as a mapping Π from 𝐺 to a set of 𝑘-tuples composed of integers Π(𝑣)={𝜋₁(𝑣),𝜋₂(𝑣),…,𝜋_𝑘(𝑣)(𝑣)}, where 𝜋_𝑖(𝑣)≤𝜋_𝑖+1(𝑣) for all 1≤𝑖<𝑘(𝑣). The simplest extension would be to introduce the case with identical sizes of 𝑘-tuples composed of nonnegative integers. A uniform  𝑘-tuple integer labeling Π_𝑐 of 𝐺 is characterized by a constant 𝑘(𝑣) for all 𝑣. The stretch of the edge 𝑣𝑤 in a Π_𝑐(𝐺) is a simple sum 𝑐_𝑣𝑤=∑^𝑘_𝑖=1𝜑(𝜋_𝑖(𝑣),𝜋_𝑖(𝑤))⋅𝜑(𝑥,𝑦) is defined as in (1). Given a uniform 𝑘-tuple integer labeling of the leaves of G the minimum sum problem is to find a labeling which minimizes the total sum of the stretches of the edges. Some 𝝅_𝐢(𝐯)=𝟎. The minimum sum problem is that of minimizing 𝐬(𝐆)=∑_{∀{𝑣𝑤}∈𝐸(𝐺)}𝑐_𝑣𝑤 over all Π_𝑐 for given 𝑘. By some modifications of the algorithms presented in this paper the minimizing 𝑘-tuple labeling can be found. This model again is a gain-loss model. More sophisticated extension must provide more realistic definition of distance between two 𝑘-tuples composed of positive integers by introducing duplication events.

Appendix

Lemma A.1 (root optimal label). Suppose that the root node is called 𝑟𝑡 and suppose that its children are called 𝑙 and 𝑟. The claim is(1)if (𝜋(𝑙)=𝜋(𝑟)=0), then 𝜋(𝑟𝑡)=0 else(2)if (𝜋(𝑙)=0), then 𝜋(𝑟𝑡)=𝜋(𝑟) else(3)if (𝜋(𝑟)=0), then 𝜋(𝑟𝑡)=𝜋(𝑙) else(4)if (𝜋(𝑙)≤𝜋(𝑟)), then 𝜋(𝑙)≤𝜋(𝑟𝑡)≤𝜋(𝑟) else(5)𝜋(𝑙)≥𝜋(𝑟𝑡)≥𝜋(𝑟).

Proof. If we consider a subtree with the node 𝑘 as a root then 𝜎_𝑘(𝑖) designates the minimal cost, given that node 𝑘 has a label 𝑖: 𝑆(𝐺)=min_𝑖𝜎_𝑁(𝑖)=min_𝑖𝜑(𝜋()(𝑙(𝜋()(𝑟𝑟𝑡,𝜋))+𝜑𝑟𝑡,𝜋))+𝜎_𝑙(𝜋(𝑙))+𝜎_𝑟(𝜋(𝑟)). (A.1) Proof of the subclaims (1)–(3) is trivial. Case (4) is 𝑆(𝐺)=min_𝑖𝜎_𝑁(𝑖)=min_𝑖|𝜋(𝑟𝑡)−𝜋(𝑙)|+|𝜋(𝑟𝑡)−𝜋(𝑟)|+𝜎_𝑙(𝜋(𝑙))+𝜎_𝑟(𝜋(𝑟)). (A.2) Let us introduce a new numbering 𝜋^′ by changing only the root label: 𝜋^′(𝑟𝑡)=𝜋(𝑙). It is easy to see that 𝑆_𝜋′(𝐺)<𝑆_𝜋(𝐺). It means that for optimal integer labeling 𝜋 the following is correct: 𝜋(𝑟𝑡)≥𝜋(𝑙). Likewise, we prove that for optimal integer labeling 𝜋(𝑟𝑡)≤𝜋(𝑟). Let us denote 𝑘=𝜋(𝑟𝑡)∶𝜋(𝑙)≤𝑘≤𝜋(𝑟): ∀𝑘(𝜋(𝑙)≤𝑘≤𝜋(𝑟))𝑆_𝜋(𝐺)=𝑘−𝜋(𝑙)+𝑆_𝜋(𝑙)+𝜋(𝑟)−𝑘+𝑆_𝜋(𝑟),q.e.d. (A.3)

Lemma A.2 (leaf parent optimal interval). Every node of the optimal integer labeling that all its children are leaf nodes has either a zero label or a label between labels of its children.
Suppose that the root node is called a and suppose that its children are called 𝑙 and 𝑟. The claim is(1)if (𝜋(𝑙)=𝜋(𝑟)=0), then 𝜋(𝑎)=0 else(2)if (𝜋(𝑙)=0), then 𝜋(𝑎)=𝜋(𝑟) else(3)if (𝜋(𝑟)=0), then 𝜋(𝑎)=𝜋(𝑙) else(4)if (𝜋(𝑙)≤𝜋(𝑟)), then 𝜋(𝑙)≤𝜋(𝑎)≤𝜋(𝑟) else(5)𝜋(𝑙)≥𝜋(𝑎)≥𝜋(𝑟).

Proof. Figure 1 illustrates this lemma. Proof of the subclaims (1)–(3) is trivial. In case of condition (4) let us prove that for optimal integer labeling 𝜋(𝑎)≥𝜋(𝑙). Suppose 𝜋(𝑎)<𝜋(𝑙). Let us denote (𝜋(𝑙)−𝜋(𝑎))=𝛿; 𝜋(𝑟)−𝜋(𝑙)=𝛾. Let us introduce a new numbering 𝜋^′ by changing only the label of the node 𝑎: 𝜋^′(𝑎)=𝜋(𝑙). It is to show that 𝑆_𝜋′(𝐺)<𝑆_𝜋(𝐺). Indeed, 𝑆_𝜋′(𝐺)=𝑆_𝜋(𝐺)−𝜋(𝑘)𝑗−𝑝−𝜋_𝑖−𝜋(𝑘)−𝑝_𝑖+1(𝑘)+𝜋−𝜋_′(𝑘)−𝜋𝑗+𝑝_𝑖−𝜋^′(𝑘)+𝑝_𝑖+1−𝜋^′(𝑘)=𝑆_𝜋(𝐺)+𝑝_𝑖−𝛿−𝜋𝑗−𝑝_𝑖𝑗−𝑝−𝜋_𝑖−𝜋(𝑘)−𝑝_𝑖+1−𝜋(𝑘)+𝑝_𝑖−𝑝_𝑖+𝑝_𝑖+1−𝑝_𝑖=𝑆_𝜋(𝐺)+𝑝_𝑖−𝛿−𝜋𝑗−𝑝_𝑖−𝜋𝑗−𝛿−𝛿+𝛾−𝛾=𝑆_𝜋(𝐺)−𝛿. (A.4)
Likewise, we prove that for optimal integer labeling 𝜋(𝑘)≤𝑝_𝑖+1. Q.e.d. 𝑝_𝑖≤𝜋(𝑘)≤𝑝_𝑖+1.

Lemma A.3 (parent optimal interval—(I)). An optimal label of a parent either is equal to zero or lies between extreme values of optimal labels of its children. If an optimal integer labeling 𝜋 provides the labels of two siblings 𝑖₁ and 𝑖₂ satisfying the conditions 𝑎₁≤𝜋(𝑖₁)≤𝑏₁ & 𝑎₂≤𝜋(𝑖₂)≤𝑏₂, then the label of their parent 𝑘 satisfies min(𝑎₁,𝑏₁,𝑎₂,𝑏₂)≤𝜋(𝑘)≤max(𝑎₁,𝑏₁,𝑎₂,𝑏₂). Proof is as for Lemma A.1.

Lemma A.4 (parent optimal interval—(II)). An optimal label of a parent in case of the empty intersection of the optimal intervals of its children lies between these intervals.
If an optimal integer labeling 𝜋 provides the labels of two siblings 𝑖₁ and 𝑖₂ satisfying the conditions (𝑎₁≤𝜋(𝑖₁)≤𝑏₁) & (𝑎₂≤𝜋(𝑖₂)≤𝑏₂) then if (𝑏₁≤𝑎₂) then 𝑏₁≤𝜋(𝑘)≤𝑎₂ else if (𝑏₂≤𝑎₁) then 𝑏₂≤𝜋(𝑘)≤𝑎₁.

Proof. (1) 𝑏₁≤𝑎₂. Let us assume 𝜋(𝑘)<𝑏₁; 𝜋(𝑘)=𝑏₁−𝛼. Then we introduce a new labeling 𝜋^′ by changing labels for three nodes: 𝜋^′(𝑘)=𝑏₁; 𝜋^′(𝑖₁)=𝑏₁; 𝜋^′(𝑖₂)=𝑎₂: 𝑆_𝜋(𝐺)−𝑆_𝜋′(𝐺)=𝜋𝑗−𝜋(𝑘)−𝜋𝑗−𝜋^′(𝑘+𝜋(𝑘)𝑖)−𝜋₁−𝜋^′(𝑘)−𝜋^′𝑖₁+𝜋(𝑘)−𝜋𝑖₂−𝜋^′(𝑘)−𝜋^′𝑖₂,𝜋𝑗−𝜋(𝑘)−𝜋𝑗−𝜋^′(𝑘=𝜋𝑗)−𝑏₁+𝛼−𝜋𝑗−𝑏₁𝜋(𝑘)𝑖=𝛼,−𝜋₁−𝜋^′(𝑘)−𝜋^′𝑖₁=𝜋(𝑘)−𝜋𝑖₁−𝜋^′(𝑘)−𝜋^′𝑖₁=𝑏₁−𝛼−𝜋𝑖₁−𝑏₁+𝑏₁=𝑏₁−𝛼−𝜋𝑖₁,𝜋(𝑘)−𝜋𝑖₂−𝜋^′(𝑘)−𝜋^′𝑖₂=𝜋𝑖₂−𝑏₁+𝛼−𝑎₂+𝑏₁,𝑆_𝜋(𝐺)−𝑆_𝜋′(𝐺)=𝛼+𝑏₁−𝛼−𝜋𝑖₁+𝜋𝑖₂+𝛼−𝑎₂=𝑏₁−𝜋𝑖₁+𝛼+𝜋𝑖₂−𝑎₂>0. (A.5) From assumption 𝜋(𝑘)<𝑏₁ follows that 𝜋 is not an optimal labeling. Similarly, we can prove that from assumption 𝜋(𝑘)>𝑎₂ follows that 𝜋 is not an optimal labeling.
(2) 𝑏₂≤𝑎₁. Similarly to (1) let us assume 𝜋(𝑘)<𝑏₂; 𝜋(𝑘)=𝑏₂−𝛼. Then we introduce a new labeling 𝜋^′ by changing labels for three nodes: 𝜋^′(𝑘)=𝑏₂; 𝜋^′(𝑖₁)=𝑎₁; 𝜋^′(𝑖₂)=𝑏₂. 𝑆_𝜋(𝐺)−𝑆_𝜋′(𝐺)>0, so we have a contradiction with the statement that 𝜋(𝐺) is an optimal labeling.

Lemma A.5 (parent optimal interval—(III)). An optimal interval of a parent is either an intersection of the optimal intervals of its children or an interval that lies between these intervals in case that their intersection is empty.
If an optimal integer labeling 𝜋 provides the labels of two siblings 𝑖₁ and 𝑖₂ satisfying the conditions 𝑎₁≤𝜋(𝑖₁)≤𝑏₁  &  𝑎₂≤𝜋(𝑖₂)≤𝑏₂, then the label of their parent 𝑘 satisfies the following condition: if ([𝑎₁,𝑏₁]∩[𝑎₂,𝑏₂])≠∅ then 𝜋(𝑘)∈[𝑎₁,𝑏₁]∩[𝑎₂,𝑏₂]) else if 𝑏₁<𝑎₂then 𝜋(𝑘)∈[𝑏₁,𝑎₂] else 𝜋(𝑘)∈[𝑏₂,𝑎₁].For example, if 𝑎₁≤𝑎₂≤𝑏₁≤𝑏₂, then Lemma A.5 states that 𝜋(𝑘) satisfies the following condition 𝑎₂≤𝜋(𝑘)≤𝑏₁. Proof is similar to proof in Lemma A.4.