--- /dev/null
+#include "muscle.h"\r
+#include "tree.h"\r
+#include <math.h>\r
+\r
+#define TRACE 0\r
+\r
+/***\r
+Sequence weights derived from a tree using Gotoh's\r
+three-way method.\r
+\r
+ Gotoh (1995) CABIOS 11(5), 543-51.\r
+\r
+Each edge e is assigned a weight w(e).\r
+\r
+Consider first a tree with three leaves A,B and C\r
+having branch lengths a, b and c, as follows.\r
+\r
+ B\r
+ |\r
+ b\r
+ |\r
+ A---a---R---c---C\r
+\r
+The internal node is denoted by R.\r
+\r
+Define:\r
+\r
+ S = (ab + ca + ab)\r
+ x = bc(a + b)(a + c)\r
+ y = a(b + c)FS\r
+\r
+Here F is a tunable normalization factor which is\r
+approximately 1.0. Then the edge weight for AR\r
+is computed as:\r
+\r
+ w(AR) = sqrt(x/y)\r
+\r
+Similar expressions for the other edges follow by\r
+symmetry.\r
+\r
+For a tree with more than three edges, the weight\r
+of an edge that ends in a leaf is computed from\r
+the three-way tree that includes the edge and\r
+its two neighbors. The weight of an internal edge\r
+is computed as the product of the weights for that\r
+edge derived from the two three-way subtrees that\r
+include that edge.\r
+\r
+For example, consider the following tree.\r
+\r
+ B\r
+ |\r
+ A--R--V--C\r
+ |\r
+ D\r
+\r
+Here, w(RV) is computed as the product of the\r
+two values for w(RV) derived from the three-way\r
+trees with leaves ABV and RCD respectively.\r
+\r
+The calculation is done using "Gotoh lengths",\r
+not the real edge lengths.\r
+\r
+The Gotoh length G of a directed edge is calculated\r
+recursively as:\r
+\r
+ G = d + LR/(L + R)\r
+\r
+where d is the length of the edge, and L and R are\r
+the Gotoh lengths of the left and right edges adjoining\r
+the terminal end of the edge. If the edge terminates on\r
+a leaf, then G=d.\r
+\r
+Pairwise sequence weights are computed as the\r
+product of edge weights on the path that connects\r
+their leaves.\r
+\r
+If the tree is split into two subtrees by deleting\r
+a given edge e, then the pairwise weights factorize.\r
+For operations on profiles formed from the two\r
+subtrees, it is possible to assign a weight to a\r
+sequence as the product of edge weights on a path\r
+from e to its leaf.\r
+***/\r
+\r
+// The xxxUnrooted functions present a rooted tree as\r
+// if it had been unrooted by deleting the root node.\r
+static unsigned GetFirstNeighborUnrooted(const Tree &tree, unsigned uNode1,\r
+ unsigned uNode2)\r
+ {\r
+ if (tree.IsRoot(uNode1) || tree.IsRoot(uNode2))\r
+ Quit("GetFirstNeighborUnrooted, should never be called with root");\r
+ if (!tree.IsEdge(uNode1, uNode2))\r
+ {\r
+ if (!tree.IsRoot(tree.GetParent(uNode1)) ||\r
+ !tree.IsRoot(tree.GetParent(uNode2)))\r
+ Quit("GetFirstNeighborUnrooted, not edge");\r
+ const unsigned uRoot = tree.GetRootNodeIndex();\r
+ return tree.GetFirstNeighbor(uNode1, uRoot);\r
+ }\r
+\r
+ unsigned uNeighbor = tree.GetFirstNeighbor(uNode1, uNode2);\r
+ if (tree.IsRoot(uNeighbor))\r
+ return tree.GetFirstNeighbor(uNeighbor, uNode1);\r
+ return uNeighbor;\r
+ }\r
+\r
+static unsigned GetSecondNeighborUnrooted(const Tree &tree, unsigned uNode1,\r
+ unsigned uNode2)\r
+ {\r
+ if (tree.IsRoot(uNode1) || tree.IsRoot(uNode2))\r
+ Quit("GetFirstNeighborUnrooted, should never be called with root");\r
+ if (!tree.IsEdge(uNode1, uNode2))\r
+ {\r
+ if (!tree.IsRoot(tree.GetParent(uNode1)) ||\r
+ !tree.IsRoot(tree.GetParent(uNode2)))\r
+ Quit("GetFirstNeighborUnrooted, not edge");\r
+ const unsigned uRoot = tree.GetRootNodeIndex();\r
+ return tree.GetSecondNeighbor(uNode1, uRoot);\r
+ }\r
+\r
+ unsigned uNeighbor = tree.GetSecondNeighbor(uNode1, uNode2);\r
+ if (tree.IsRoot(uNeighbor))\r
+ return tree.GetFirstNeighbor(uNeighbor, uNode1);\r
+ return uNeighbor;\r
+ }\r
+\r
+static unsigned GetNeighborUnrooted(const Tree &tree, unsigned uNode1,\r
+ unsigned uSub)\r
+ {\r
+ unsigned uNeighbor = tree.GetNeighbor(uNode1, uSub);\r
+ if (tree.IsRoot(uNeighbor))\r
+ return tree.GetFirstNeighbor(uNeighbor, uNode1);\r
+ return uNeighbor;\r
+ }\r
+\r
+static unsigned GetNeighborSubscriptUnrooted(const Tree &tree, unsigned uNode1,\r
+ unsigned uNode2)\r
+ {\r
+ if (tree.IsEdge(uNode1, uNode2))\r
+ return tree.GetNeighborSubscript(uNode1, uNode2);\r
+ if (!tree.IsRoot(tree.GetParent(uNode1)) ||\r
+ !tree.IsRoot(tree.GetParent(uNode2)))\r
+ Quit("GetNeighborSubscriptUnrooted, not edge");\r
+ for (unsigned uSub = 0; uSub < 3; ++uSub)\r
+ if (GetNeighborUnrooted(tree, uNode1, uSub) == uNode2)\r
+ return uSub;\r
+ Quit("GetNeighborSubscriptUnrooted, not a neighbor");\r
+ return NULL_NEIGHBOR;\r
+ }\r
+\r
+static double GetEdgeLengthUnrooted(const Tree &tree, unsigned uNode1,\r
+ unsigned uNode2)\r
+ {\r
+ if (tree.IsRoot(uNode1) || tree.IsRoot(uNode2))\r
+ Quit("GetEdgeLengthUnrooted, should never be called with root");\r
+ if (!tree.IsEdge(uNode1, uNode2))\r
+ {\r
+ if (!tree.IsRoot(tree.GetParent(uNode1)) ||\r
+ !tree.IsRoot(tree.GetParent(uNode2)))\r
+ Quit("GetEdgeLengthUnrooted, not edge");\r
+\r
+ const unsigned uRoot = tree.GetRootNodeIndex();\r
+ return tree.GetEdgeLength(uNode1, uRoot) +\r
+ tree.GetEdgeLength(uNode2, uRoot);\r
+ }\r
+ return tree.GetEdgeLength(uNode1, uNode2);\r
+ }\r
+\r
+double GetGotohLength(const Tree &tree, unsigned R, unsigned A)\r
+ {\r
+ double dThis = GetEdgeLengthUnrooted(tree, R, A);\r
+\r
+// Enforce non-negative edge lengths\r
+ if (dThis < 0)\r
+ dThis = 0;\r
+\r
+ if (tree.IsLeaf(A))\r
+ return dThis;\r
+\r
+ const unsigned uFirst = GetFirstNeighborUnrooted(tree, A, R);\r
+ const unsigned uSecond = GetSecondNeighborUnrooted(tree, A, R);\r
+ const double dFirst = GetGotohLength(tree, A, uFirst);\r
+ const double dSecond = GetGotohLength(tree, A, uSecond);\r
+ const double dSum = dFirst + dSecond;\r
+ const double dThird = dSum == 0 ? 0 : (dFirst*dSecond)/dSum;\r
+ return dThis + dThird;\r
+ }\r
+\r
+// Return weight of edge A-R in three-way subtree that has\r
+// leaves A,B,C and internal node R.\r
+static double GotohWeightThreeWay(const Tree &tree, unsigned A,\r
+ unsigned B, unsigned C, unsigned R)\r
+ {\r
+ const double F = 1.0;\r
+\r
+ if (tree.IsLeaf(R))\r
+ Quit("GotohThreeWay: R must be internal node");\r
+\r
+ double a = GetGotohLength(tree, R, A);\r
+ double b = GetGotohLength(tree, R, B);\r
+ double c = GetGotohLength(tree, R, C);\r
+\r
+ double S = b*c + c*a + a*b;\r
+ double x = b*c*(a + b)*(a + c);\r
+ double y = a*(b + c)*F*S;\r
+\r
+// y is zero iff all three branch lengths are zero.\r
+ if (y < 0.001)\r
+ return 1.0;\r
+ return sqrt(x/y);\r
+ }\r
+\r
+static double GotohWeightEdge(const Tree &tree, unsigned uNodeIndex1,\r
+ unsigned uNodeIndex2)\r
+ {\r
+ double w1 = 1.0;\r
+ double w2 = 1.0;\r
+ if (!tree.IsLeaf(uNodeIndex1))\r
+ {\r
+ unsigned R = uNodeIndex1;\r
+ unsigned A = uNodeIndex2;\r
+ unsigned B = GetFirstNeighborUnrooted(tree, R, A);\r
+ unsigned C = GetSecondNeighborUnrooted(tree, R, A);\r
+ w1 = GotohWeightThreeWay(tree, A, B, C, R);\r
+ }\r
+ if (!tree.IsLeaf(uNodeIndex2))\r
+ {\r
+ unsigned R = uNodeIndex2;\r
+ unsigned A = uNodeIndex1;\r
+ unsigned B = GetFirstNeighborUnrooted(tree, R, A);\r
+ unsigned C = GetSecondNeighborUnrooted(tree, R, A);\r
+ w2 = GotohWeightThreeWay(tree, A, B, C, R);\r
+ }\r
+ return w1*w2;\r
+ }\r
+\r
+void CalcThreeWayEdgeWeights(const Tree &tree, WEIGHT **EdgeWeights)\r
+ {\r
+ const unsigned uNodeCount = tree.GetNodeCount();\r
+ for (unsigned uNodeIndex1 = 0; uNodeIndex1 < uNodeCount; ++uNodeIndex1)\r
+ {\r
+ if (tree.IsRoot(uNodeIndex1))\r
+ continue;\r
+ for (unsigned uSub1 = 0; uSub1 < 3; ++uSub1)\r
+ {\r
+ const unsigned uNodeIndex2 = GetNeighborUnrooted(tree, uNodeIndex1, uSub1);\r
+ if (NULL_NEIGHBOR == uNodeIndex2)\r
+ continue;\r
+\r
+ // Avoid computing same edge twice in reversed order\r
+ if (uNodeIndex2 < uNodeIndex1)\r
+ continue;\r
+\r
+ const WEIGHT w = (WEIGHT) GotohWeightEdge(tree, uNodeIndex1, uNodeIndex2);\r
+ const unsigned uSub2 = GetNeighborSubscriptUnrooted(tree, uNodeIndex2, uNodeIndex1);\r
+#if DEBUG\r
+ {\r
+ assert(uNodeIndex2 == GetNeighborUnrooted(tree, uNodeIndex1, uSub1));\r
+ assert(uNodeIndex1 == GetNeighborUnrooted(tree, uNodeIndex2, uSub2));\r
+ const WEIGHT wRev = (WEIGHT) GotohWeightEdge(tree, uNodeIndex2, uNodeIndex1);\r
+ if (!BTEq(w, wRev))\r
+ Quit("CalcThreeWayWeights: rev check failed %g %g",\r
+ w, wRev);\r
+ }\r
+#endif\r
+ EdgeWeights[uNodeIndex1][uSub1] = w;\r
+ EdgeWeights[uNodeIndex2][uSub2] = w;\r
+ }\r
+ }\r
+ }\r
+\r
+static void SetSeqWeights(const Tree &tree, unsigned uNode1, unsigned uNode2,\r
+ double dPathWeight, WEIGHT *Weights)\r
+ {\r
+ if (tree.IsRoot(uNode1) || tree.IsRoot(uNode2))\r
+ Quit("SetSeqWeights, should never be called with root");\r
+\r
+ const double dThisLength = GetEdgeLengthUnrooted(tree, uNode1, uNode2);\r
+ if (tree.IsLeaf(uNode2))\r
+ {\r
+ const unsigned Id = tree.GetLeafId(uNode2);\r
+ Weights[Id] = (WEIGHT) (dPathWeight + dThisLength);\r
+ return;\r
+ }\r
+ const unsigned uFirst = GetFirstNeighborUnrooted(tree, uNode2, uNode1);\r
+ const unsigned uSecond = GetSecondNeighborUnrooted(tree, uNode2, uNode1);\r
+ dPathWeight *= dThisLength;\r
+ SetSeqWeights(tree, uNode2, uFirst, dPathWeight, Weights);\r
+ SetSeqWeights(tree, uNode2, uSecond, dPathWeight, Weights);\r
+ }\r
+\r
+void CalcThreeWayWeights(const Tree &tree, unsigned uNode1, unsigned uNode2,\r
+ WEIGHT *Weights)\r
+ {\r
+#if TRACE\r
+ Log("CalcThreeWayEdgeWeights\n");\r
+ tree.LogMe();\r
+#endif\r
+\r
+ if (tree.IsRoot(uNode1))\r
+ uNode1 = tree.GetFirstNeighbor(uNode1, uNode2);\r
+ else if (tree.IsRoot(uNode2))\r
+ uNode2 = tree.GetFirstNeighbor(uNode2, uNode1);\r
+ const unsigned uNodeCount = tree.GetNodeCount();\r
+ WEIGHT **EdgeWeights = new WEIGHT *[uNodeCount];\r
+ for (unsigned uNodeIndex = 0; uNodeIndex < uNodeCount; ++uNodeIndex)\r
+ EdgeWeights[uNodeIndex] = new WEIGHT[3];\r
+\r
+ CalcThreeWayEdgeWeights(tree, EdgeWeights);\r
+\r
+#if TRACE\r
+ {\r
+ Log("Node1 Node2 Length Gotoh EdgeWt\n");\r
+ Log("----- ----- ------ ------ ------\n");\r
+ for (unsigned uNodeIndex1 = 0; uNodeIndex1 < uNodeCount; ++uNodeIndex1)\r
+ {\r
+ if (tree.IsRoot(uNodeIndex1))\r
+ continue;\r
+ for (unsigned uSub1 = 0; uSub1 < 3; ++uSub1)\r
+ {\r
+ const unsigned uNodeIndex2 = GetNeighborUnrooted(tree, uNodeIndex1, uSub1);\r
+ if (NULL_NEIGHBOR == uNodeIndex2)\r
+ continue;\r
+ if (uNodeIndex2 < uNodeIndex1)\r
+ continue;\r
+ const WEIGHT ew = EdgeWeights[uNodeIndex1][uSub1];\r
+ const double d = GetEdgeLengthUnrooted(tree, uNodeIndex1, uNodeIndex2);\r
+ const double g = GetGotohLength(tree, uNodeIndex1, uNodeIndex2);\r
+ Log("%5u %5u %6.3f %6.3f %6.3f\n", uNodeIndex1, uNodeIndex2, d, g, ew);\r
+ }\r
+ }\r
+ }\r
+#endif\r
+\r
+ SetSeqWeights(tree, uNode1, uNode2, 0.0, Weights);\r
+ SetSeqWeights(tree, uNode2, uNode1, 0.0, Weights);\r
+\r
+ for (unsigned uNodeIndex = 0; uNodeIndex < uNodeCount; ++uNodeIndex)\r
+ delete[] EdgeWeights[uNodeIndex];\r
+ delete[] EdgeWeights;\r
+ }\r
git://source.jalview.org / jabaws.git / blobdiff