--- /dev/null
+/***\r
+Code for implementing HMMer's "BLOSUM weighting" algorithm.\r
+\r
+The algorithm was deduced by reverse-engineering the HMMer code.\r
+\r
+The HMMer documentation refers to BLOSUM weighting as "Henikoff\r
+simple filter weighting"\r
+\r
+The name BLOSUM implied to me that HMMer would be using a\r
+substitution probability matrix to compute distances, but this\r
+turned out not to be the case.\r
+\r
+It is notable, not to say puzzling, that the HMMer BLOSUM weighting\r
+algorithm is guaranteed to produce an integral NIC (number-of-indepdent-\r
+counts, also known as effective sequence count). Presumably Eddy must\r
+have known this, though he doesn't comment on it and he computes & stores\r
+the value in a float.\r
+\r
+Here's the algorithm:\r
+\r
+Distances between two sequences are based on the average of a simple \r
+binary equal (one) / not equal (zero) at each position. The only thing\r
+that has anything to do with BLOSUM in this calculation is an obscure\r
+(to me) constant value of 0.62. The sequences are clustered using this\r
+distance. If the pairwise identity (fraction of identical positions)\r
+is less than 0.62, they get assigned to disjoint clusters, the final\r
+number of disjoint clusters is the NIC. This makes some intuitive sense:\r
+I would interpret this by saying that if a set of sequences are close\r
+enough they count as one sequence. The weight for each sequence within a\r
+disjoint cluster is then determined to be 1 / (clustersize), from which it\r
+follows that the sum of all weights is equal to the number of disjoint\r
+clusters and is thus guaranteed to be an integer value. It is exactly this\r
+sum that HMMer uses for the NIC, by default.\r
+\r
+The individual BLOSUM sequence weights are not used for anything else in\r
+HMMer, unless you specify that BLOSUM weighting should override the default\r
+GSC weighting. GSC weighting uses a different clustering algorithm to\r
+determine relative weights. The BLOSUM NIC is then distributed over the\r
+GSC tree according to those relative weights.\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::SetBLOSUMSubtreeWeight(const ClusterNode *ptrNode, double dWeight) 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
+ WEIGHT w = DoubleToWeight(dWeight);\r
+ m_Weights[uIndex] = w;\r
+ return;\r
+ }\r
+\r
+// Otherwise, recursively set subtrees\r
+ SetBLOSUMSubtreeWeight(ptrLeft, dWeight);\r
+ SetBLOSUMSubtreeWeight(ptrRight, dWeight);\r
+ }\r
+\r
+// Traverse a subtree looking for clusters where all\r
+// the leaves are sufficiently similar that they\r
+// should be weighted as a group, i.e. given a weight\r
+// of 1/N where N is the cluster size. The idea is\r
+// to avoid sample bias where we have closely related\r
+// sequences in the input alignment.\r
+// The weight at a node is the distance between\r
+// the two closest sequences in the left and right\r
+// subtrees under that node. "Sufficiently similar"\r
+// is defined as being where that minimum distance\r
+// is less than the dMinDist threshhold. I don't know\r
+// why the clustering is done using a minimum rather\r
+// than a maximum or average, either of which would\r
+// seem more natural to me.\r
+// Return value is number of groups under this node.\r
+// A "group" is the cluster found under a node with a\r
+// weight less than the minimum.\r
+unsigned MSA::SetBLOSUMNodeWeight(const ClusterNode *ptrNode, double dMinDist) const\r
+ {\r
+ if (0 == ptrNode)\r
+ return 0;\r
+\r
+ if (ptrNode->GetWeight() < dMinDist)\r
+ {\r
+ unsigned uClusterSize = ptrNode->GetClusterSize(); \r
+ assert(uClusterSize > 0);\r
+ double dWeight = 1.0 / uClusterSize;\r
+ SetBLOSUMSubtreeWeight(ptrNode, dWeight);\r
+ return 1;\r
+ }\r
+\r
+ const ClusterNode *ptrLeft = ptrNode->GetLeft();\r
+ const ClusterNode *ptrRight = ptrNode->GetRight();\r
+\r
+ unsigned uLeftGroupCount = SetBLOSUMNodeWeight(ptrLeft, dMinDist);\r
+ unsigned uRightGroupCount = SetBLOSUMNodeWeight(ptrRight, dMinDist);\r
+\r
+ return uLeftGroupCount + uRightGroupCount;\r
+ }\r
+\r
+// Return value is the group count, i.e. the effective number\r
+// of distinctly different sequences.\r
+unsigned MSA::CalcBLOSUMWeights(ClusterTree &BlosumCluster) const\r
+ {\r
+// Build distance matrix\r
+ DistFunc DF;\r
+ unsigned uSeqCount = GetSeqCount();\r
+ DF.SetCount(uSeqCount);\r
+ for (unsigned i = 0; i < uSeqCount; ++i)\r
+ for (unsigned j = i+1; j < uSeqCount; ++j)\r
+ {\r
+ double dDist = GetPctIdentityPair(i, j);\r
+ assert(dDist >= 0.0 && dDist <= 1.0);\r
+ DF.SetDist(i, j, (float) (1.0 - dDist));\r
+ }\r
+\r
+// Cluster based on the distance function\r
+ BlosumCluster.Create(DF);\r
+\r
+// Return value is HMMer's "effective sequence count".\r
+ return SetBLOSUMNodeWeight(BlosumCluster.GetRoot(), 1.0 - BLOSUM_DIST);\r
+ }\r
git://source.jalview.org / jabaws.git / blobdiff