binaries/src/muscle/setblosumweights.cpp

   1 /***\r
   2 Code for implementing HMMer's "BLOSUM weighting" algorithm.\r
   3 \r
   4 The algorithm was deduced by reverse-engineering the HMMer code.\r
   5 \r
   6 The HMMer documentation refers to BLOSUM weighting as "Henikoff\r
   7 simple filter weighting"\r
   8 \r
   9 The name BLOSUM implied to me that HMMer would be using a\r
  10 substitution probability matrix to compute distances, but this\r
  11 turned out not to be the case.\r
  12 \r
  13 It is notable, not to say puzzling, that the HMMer BLOSUM weighting\r
  14 algorithm  is guaranteed to produce an integral NIC (number-of-indepdent-\r
  15 counts, also known as effective sequence count). Presumably Eddy must\r
  16 have known this, though he doesn't comment on it and he computes & stores\r
  17 the value in a float.\r
  18 \r
  19 Here's the algorithm:\r
  20 \r
  21 Distances between two sequences are based on the average of a simple \r
  22 binary equal (one) / not equal (zero) at each position. The only thing\r
  23 that has  anything to do with BLOSUM in this calculation is an obscure\r
  24 (to me) constant  value of 0.62. The sequences are clustered using this\r
  25 distance. If the pairwise identity (fraction of  identical positions)\r
  26 is less than 0.62, they get assigned to disjoint clusters, the final\r
  27 number of disjoint clusters is the NIC. This makes some intuitive sense:\r
  28 I would interpret this by saying that if a set of sequences are close\r
  29 enough they count as one sequence. The weight for each sequence within a\r
  30 disjoint cluster is then determined to be 1 / (clustersize), from which it\r
  31 follows that the sum of all weights is equal to the number of disjoint\r
  32 clusters and is thus guaranteed to be an integer value. It is exactly this\r
  33 sum that HMMer uses for the NIC, by default.\r
  34 \r
  35 The individual BLOSUM sequence weights are not used for anything else in\r
  36 HMMer, unless you specify that BLOSUM weighting should override the default\r
  37 GSC  weighting. GSC weighting uses a different clustering algorithm to\r
  38 determine  relative weights. The BLOSUM NIC is then distributed over the\r
  39 GSC tree according to those relative weights.\r
  40 ***/\r
  41 \r
  42 #include "muscle.h"\r
  43 #include "msa.h"\r
  44 #include "cluster.h"\r
  45 #include "distfunc.h"\r
  46 \r
  47 // Set weights of all sequences in the subtree under given node.\r
  48 void MSA::SetBLOSUMSubtreeWeight(const ClusterNode *ptrNode, double dWeight) const\r
  49         {\r
  50         if (0 == ptrNode)\r
  51                 return;\r
  52 \r
  53         const ClusterNode *ptrRight = ptrNode->GetRight();\r
  54         const ClusterNode *ptrLeft = ptrNode->GetLeft();\r
  55 \r
  56 // If leaf, set weight\r
  57         if (0 == ptrRight && 0 == ptrLeft)\r
  58                 {\r
  59                 unsigned uIndex = ptrNode->GetIndex();\r
  60                 WEIGHT w = DoubleToWeight(dWeight);\r
  61                 m_Weights[uIndex] = w;\r
  62                 return;\r
  63                 }\r
  64 \r
  65 // Otherwise, recursively set subtrees\r
  66         SetBLOSUMSubtreeWeight(ptrLeft, dWeight);\r
  67         SetBLOSUMSubtreeWeight(ptrRight, dWeight);\r
  68         }\r
  69 \r
  70 // Traverse a subtree looking for clusters where all\r
  71 // the leaves are sufficiently similar that they\r
  72 // should be weighted as a group, i.e. given a weight\r
  73 // of 1/N where N is the cluster size. The idea is\r
  74 // to avoid sample bias where we have closely related\r
  75 // sequences in the input alignment.\r
  76 // The weight at a node is the distance between\r
  77 // the two closest sequences in the left and right\r
  78 // subtrees under that node. "Sufficiently similar"\r
  79 // is defined as being where that minimum distance\r
  80 // is less than the dMinDist threshhold. I don't know\r
  81 // why the clustering is done using a minimum rather\r
  82 // than a maximum or average, either of which would\r
  83 // seem more natural to me.\r
  84 // Return value is number of groups under this node.\r
  85 // A "group" is the cluster found under a node with a\r
  86 // weight less than the minimum.\r
  87 unsigned MSA::SetBLOSUMNodeWeight(const ClusterNode *ptrNode, double dMinDist) const\r
  88         {\r
  89         if (0 == ptrNode)\r
  90                 return 0;\r
  91 \r
  92         if (ptrNode->GetWeight() < dMinDist)\r
  93                 {\r
  94                 unsigned uClusterSize = ptrNode->GetClusterSize(); \r
  95                 assert(uClusterSize > 0);\r
  96                 double dWeight = 1.0 / uClusterSize;\r
  97                 SetBLOSUMSubtreeWeight(ptrNode, dWeight);\r
  98                 return 1;\r
  99                 }\r
 100 \r
 101         const ClusterNode *ptrLeft = ptrNode->GetLeft();\r
 102         const ClusterNode *ptrRight = ptrNode->GetRight();\r
 103 \r
 104         unsigned uLeftGroupCount = SetBLOSUMNodeWeight(ptrLeft, dMinDist);\r
 105         unsigned uRightGroupCount = SetBLOSUMNodeWeight(ptrRight, dMinDist);\r
 106 \r
 107         return uLeftGroupCount + uRightGroupCount;\r
 108         }\r
 109 \r
 110 // Return value is the group count, i.e. the effective number\r
 111 // of distinctly different sequences.\r
 112 unsigned MSA::CalcBLOSUMWeights(ClusterTree &BlosumCluster) const\r
 113         {\r
 114 // Build distance matrix\r
 115         DistFunc DF;\r
 116         unsigned uSeqCount = GetSeqCount();\r
 117         DF.SetCount(uSeqCount);\r
 118         for (unsigned i = 0; i < uSeqCount; ++i)\r
 119                 for (unsigned j = i+1; j < uSeqCount; ++j)\r
 120                         {\r
 121                         double dDist = GetPctIdentityPair(i, j);\r
 122                         assert(dDist >= 0.0 && dDist <= 1.0);\r
 123                         DF.SetDist(i, j, (float) (1.0 - dDist));\r
 124                         }\r
 125 \r
 126 // Cluster based on the distance function\r
 127         BlosumCluster.Create(DF);\r
 128 \r
 129 // Return value is HMMer's "effective sequence count".\r
 130         return SetBLOSUMNodeWeight(BlosumCluster.GetRoot(), 1.0 - BLOSUM_DIST);\r
 131         }\r