+++ /dev/null
-/***\r
-Gerstein/Sonnhammer/Chothia ad hoc sequence weighting.\r
-The algorithm was deduced by reverse-engineering the\r
-HMMer code.\r
-\r
-I used an alternative representation that I prefer over\r
-HMMer's. The HMMer code is full of tree manipulations\r
-that do something to the left child and then the equivalent\r
-thing to the right child. It was clear that there must be\r
-a re-formulation that does everything once for each node,\r
-which would reduce the number of operations expressed\r
-in the code by a factor of two. This gives a more elegant\r
-and less error-prone way to code it.\r
-\r
-These notes explain the correspondence between my design\r
-and Eddy's.\r
-\r
-HMMer stores a data structure phylo_s for each non-leaf\r
-node in the cluster tree. This structure contains the\r
-following fields:\r
-\r
- diff Weight of the node\r
- lblen Left branch length\r
- rblen Right branch length\r
-\r
-The lblen and rblen branch lengths are calculated as:\r
-\r
- this.lblen = this.diff - left.diff\r
- this.rblen = this.diff - right.diff\r
-\r
-My code stores one ClusterNode data structure per node\r
-in the cluster tree, including leaves. I store only the\r
-weight. I can recover the HMMer branch length fields\r
-in a trivial O(1) calculation as follows:\r
-\r
- lblen = Node.GetWeight() - Node.GetLeft()->GetWeight()\r
- rblen = Node.GetWeight() - Node.GetRight()->GetWeight()\r
-\r
-For the GSC weights calculation, HMMer constructs the\r
-following vectors, which have entries for all nodes,\r
-including leaves:\r
-\r
- lwt Left weight\r
- rwt Right weight\r
-\r
-The "left weight" is calculated as the sum of the weights in\r
-all the nodes reachable through the left branch, including\r
-the node itself. (This is not immediately obvious from the\r
-code, which does the calculation using branch lengths rather\r
-than weights, but this is an equivalent, and to my mind clearer,\r
-statement of what they are). Similarly, the "right weight" is\r
-the sum of all weights reachable via the right branch. I define\r
-the "cluster weight" to be the summed weight of all nodes in the\r
-subtree under the node, including the node itself. I provide\r
-a function Node.GetClusterWeight() which calculates the cluster\r
-weight using a O(ln N) recursion through the tree. The lwt and\r
-rwt values can be recovered as follows:\r
-\r
- lwt = Node.GetLeft()->GetClusterWeight()\r
- + Node.GetWeight()\r
-\r
- lwt = Node.GetLeft()->GetClusterWeight()\r
- + Node.GetWeight()\r
-\r
-HMMer calculates a further vector fwt as follows.\r
-\r
- this.fwt = parent.fwt * parent.lwt / (parent.lwt + parent.rwt)\r
-\r
-This applies to nodes reached via a left branch, for nodes reached\r
-via a right branch:\r
-\r
- this.fwt = parent.fwt * parent.rwt / (parent.lwt + parent.rwt)\r
-\r
-The values of fwt at the leaf nodes are the final GSC weights.\r
-We derive the various terms using our equivalents.\r
-\r
- parent.lwt = Parent.GetLeft()->GetClusterWeight()\r
- + Parent.GetWeight()\r
-\r
- parent.rwt = Parent.GetRight()->GetClusterWeight()\r
- + Parent.GetWeight()\r
-\r
- parent.lwt + parent.rwt =\r
- { Parent.GetLeft()->GetClusterWeight()\r
- + Parent.GetRight()->GetClusterWeight()\r
- + Parent.GetWeight() }\r
- + Parent.GetWeight()\r
-\r
-We recognize the term {...} as the cluster weight of the\r
-parent, so\r
-\r
- parent.lwt + parent.rwt\r
- = Parent.GetClusterWeight()\r
- + Parent.GetWeight()\r
-\r
-As you would expect, repeating this exercise for parent.rwt gives\r
-exactly the same expression.\r
-\r
-The GSC weights (fwt) are stored in the Weight2 field of the cluster\r
-tree, the Weight field stores the original (BLOSUM) weights used\r
-as input to this algorithm.\r
-***/\r
-\r
-#include "muscle.h"\r
-#include "msa.h"\r
-#include "cluster.h"\r
-#include "distfunc.h"\r
-\r
-// Set weights of all sequences in the subtree under given node.\r
-void MSA::SetSubtreeWeight2(const ClusterNode *ptrNode) const\r
- {\r
- if (0 == ptrNode)\r
- return;\r
-\r
- const ClusterNode *ptrRight = ptrNode->GetRight();\r
- const ClusterNode *ptrLeft = ptrNode->GetLeft();\r
-\r
-// If leaf, set weight\r
- if (0 == ptrRight && 0 == ptrLeft)\r
- {\r
- unsigned uIndex = ptrNode->GetIndex();\r
- double dWeight = ptrNode->GetWeight2();\r
- WEIGHT w = DoubleToWeight(dWeight);\r
- m_Weights[uIndex] = w;\r
- return;\r
- }\r
-\r
-// Otherwise, recursively set subtrees\r
- SetSubtreeWeight2(ptrLeft);\r
- SetSubtreeWeight2(ptrRight);\r
- }\r
-\r
-void MSA::SetSubtreeGSCWeight(ClusterNode *ptrNode) const\r
- {\r
- if (0 == ptrNode)\r
- return;\r
-\r
- ClusterNode *ptrParent = ptrNode->GetParent();\r
- double dParentWeight2 = ptrParent->GetWeight2();\r
- double dParentClusterWeight = ptrParent->GetClusterWeight();\r
- if (0.0 == dParentClusterWeight)\r
- {\r
- double dThisClusterSize = ptrNode->GetClusterSize();\r
- double dParentClusterSize = ptrParent->GetClusterSize();\r
- double dWeight2 =\r
- dParentWeight2*dThisClusterSize/dParentClusterSize;\r
- ptrNode->SetWeight2(dWeight2);\r
- }\r
- else\r
- {\r
- // Could cache cluster weights for better performance.\r
- // We calculate cluster weight of each node twice, so this\r
- // would give x2 improvement.\r
- // As weighting is not very expensive, we don't care.\r
- double dThisClusterWeight = ptrNode->GetClusterWeight();\r
- double dParentWeight = ptrParent->GetWeight();\r
-\r
- double dNum = dThisClusterWeight + dParentWeight;\r
- double dDenom = dParentClusterWeight + dParentWeight;\r
- double dWeight2 = dParentWeight2*(dNum/dDenom);\r
-\r
- ptrNode->SetWeight2(dWeight2);\r
- }\r
-\r
- SetSubtreeGSCWeight(ptrNode->GetLeft());\r
- SetSubtreeGSCWeight(ptrNode->GetRight());\r
- }\r
-\r
-void MSA::SetGSCWeights() const\r
- {\r
- ClusterTree CT;\r
- CalcBLOSUMWeights(CT);\r
-\r
-// Calculate weights and store in tree.\r
- ClusterNode *ptrRoot = CT.GetRoot();\r
- ptrRoot->SetWeight2(1.0);\r
- SetSubtreeGSCWeight(ptrRoot->GetLeft());\r
- SetSubtreeGSCWeight(ptrRoot->GetRight());\r
-\r
-// Copy weights from tree to MSA.\r
- SetSubtreeWeight2(ptrRoot);\r
- }\r
- \r
-void MSA::ListWeights() const\r
- {\r
- const unsigned uSeqCount = GetSeqCount();\r
- Log("Weights:\n");\r
- WEIGHT wTotal = 0;\r
- for (unsigned n = 0; n < uSeqCount; ++n)\r
- {\r
- wTotal += GetSeqWeight(n);\r
- Log("%6.3f %s\n", GetSeqWeight(n), GetSeqName(n));\r
- }\r
- Log("Total weights = %6.3f, should be 1.0\n", wTotal);\r
- }\r
git://source.jalview.org / jabaws.git / blobdiff