binaries/src/muscle/realigndiffse.cpp

   1 #include "muscle.h"\r
   2 #include "msa.h"\r
   3 #include "tree.h"\r
   4 #include "profile.h"\r
   5 #include "pwpath.h"\r
   6 #include "seqvect.h"\r
   7 #include "estring.h"\r
   8 \r
   9 #define TRACE           0\r
  10 \r
  11 void DeleteProgNode(ProgNode &Node)\r
  12         {\r
  13         delete[] Node.m_Prof;\r
  14         delete[] Node.m_EstringL;\r
  15         delete[] Node.m_EstringR;\r
  16 \r
  17         Node.m_Prof = 0;\r
  18         Node.m_EstringL = 0;\r
  19         Node.m_EstringR = 0;\r
  20         }\r
  21 \r
  22 static void MakeNode(ProgNode &OldNode, ProgNode &NewNode, bool bSwapLR)\r
  23         {\r
  24         if (bSwapLR)\r
  25                 {\r
  26                 NewNode.m_EstringL = OldNode.m_EstringR;\r
  27                 NewNode.m_EstringR = OldNode.m_EstringL;\r
  28                 }\r
  29         else\r
  30                 {\r
  31                 NewNode.m_EstringL = OldNode.m_EstringL;\r
  32                 NewNode.m_EstringR = OldNode.m_EstringR;\r
  33                 }\r
  34         NewNode.m_Prof = OldNode.m_Prof;\r
  35         NewNode.m_uLength = OldNode.m_uLength;\r
  36         NewNode.m_Weight = OldNode.m_Weight;\r
  37 \r
  38         OldNode.m_Prof = 0;\r
  39         OldNode.m_EstringL = 0;\r
  40         OldNode.m_EstringR = 0;\r
  41         }\r
  42 \r
  43 void RealignDiffsE(const MSA &msaIn, const SeqVect &v,\r
  44   const Tree &NewTree, const Tree &OldTree, \r
  45   const unsigned uNewNodeIndexToOldNodeIndex[],\r
  46   MSA &msaOut, ProgNode *OldProgNodes)\r
  47         {\r
  48         assert(OldProgNodes != 0);\r
  49 \r
  50         const unsigned uNodeCount = NewTree.GetNodeCount();\r
  51         if (uNodeCount%2 == 0)\r
  52                 Quit("RealignDiffs: Expected odd number of nodes");\r
  53 \r
  54         const unsigned uMergeCount = (uNodeCount - 1)/2;\r
  55         ProgNode *NewProgNodes = new ProgNode[uNodeCount];\r
  56 \r
  57         for (unsigned uNewNodeIndex = 0; uNewNodeIndex < uNodeCount; ++uNewNodeIndex)\r
  58                 {\r
  59                 if (NODE_CHANGED == uNewNodeIndexToOldNodeIndex[uNewNodeIndex])\r
  60                         continue;\r
  61 \r
  62                 unsigned uOldNodeIndex = uNewNodeIndexToOldNodeIndex[uNewNodeIndex];\r
  63                 assert(uNewNodeIndex < uNodeCount);\r
  64                 assert(uOldNodeIndex < uNodeCount);\r
  65 \r
  66                 ProgNode &NewNode = NewProgNodes[uNewNodeIndex];\r
  67                 ProgNode &OldNode = OldProgNodes[uOldNodeIndex];\r
  68                 bool bSwapLR = false;\r
  69                 if (!NewTree.IsLeaf(uNewNodeIndex))\r
  70                         {\r
  71                         unsigned uNewLeft = NewTree.GetLeft(uNewNodeIndex);\r
  72                         unsigned uNewRight = NewTree.GetRight(uNewNodeIndex);\r
  73                         unsigned uOld = uNewNodeIndexToOldNodeIndex[uNewNodeIndex];\r
  74                         unsigned uOldLeft = OldTree.GetLeft(uOld);\r
  75                         unsigned uOldRight = OldTree.GetRight(uOld);\r
  76                         assert(uOldLeft < uNodeCount && uOldRight < uNodeCount);\r
  77                         if (uOldLeft != uNewNodeIndexToOldNodeIndex[uNewLeft])\r
  78                                 {\r
  79                                 assert(uOldLeft == uNewNodeIndexToOldNodeIndex[uNewRight]);\r
  80                                 bSwapLR = true;\r
  81                                 }\r
  82                         }\r
  83                 MakeNode(OldNode, NewNode, bSwapLR);\r
  84 #if     TRACE\r
  85                 Log("MakeNode old=%u new=%u swap=%d length=%u weight=%.3g\n",\r
  86                   uOldNodeIndex, uNewNodeIndex, bSwapLR, NewNode.m_uLength, NewNode.m_Weight);\r
  87 #endif\r
  88                 }\r
  89 \r
  90         unsigned uJoin = 0;\r
  91         SetProgressDesc("Refine tree");\r
  92         for (unsigned uNewNodeIndex = NewTree.FirstDepthFirstNode();\r
  93           NULL_NEIGHBOR != uNewNodeIndex;\r
  94           uNewNodeIndex = NewTree.NextDepthFirstNode(uNewNodeIndex))\r
  95                 {\r
  96                 if (NODE_CHANGED != uNewNodeIndexToOldNodeIndex[uNewNodeIndex])\r
  97                         continue;\r
  98 \r
  99                 Progress(uJoin, uMergeCount - 1);\r
 100                 ++uJoin;\r
 101 \r
 102                 const unsigned uMergeNodeIndex = uNewNodeIndex;\r
 103                 ProgNode &Parent = NewProgNodes[uMergeNodeIndex];\r
 104 \r
 105                 const unsigned uLeft = NewTree.GetLeft(uNewNodeIndex);\r
 106                 const unsigned uRight = NewTree.GetRight(uNewNodeIndex);\r
 107 \r
 108                 ProgNode &Node1 = NewProgNodes[uLeft];\r
 109                 ProgNode &Node2 = NewProgNodes[uRight];\r
 110 \r
 111                 AlignTwoProfs(\r
 112                         Node1.m_Prof, Node1.m_uLength, Node1.m_Weight,\r
 113                         Node2.m_Prof, Node2.m_uLength, Node2.m_Weight,\r
 114                         Parent.m_Path,\r
 115                         &Parent.m_Prof, &Parent.m_uLength);\r
 116                 PathToEstrings(Parent.m_Path, &Parent.m_EstringL, &Parent.m_EstringR);\r
 117 \r
 118                 Parent.m_Weight = Node1.m_Weight + Node2.m_Weight;\r
 119 \r
 120                 delete[] Node1.m_Prof;\r
 121                 delete[] Node2.m_Prof;\r
 122 \r
 123                 Node1.m_Prof = 0;\r
 124                 Node2.m_Prof = 0;\r
 125                 }\r
 126 \r
 127         ProgressStepsDone();\r
 128 \r
 129         if (g_bBrenner)\r
 130                 MakeRootMSABrenner((SeqVect &) v, NewTree, NewProgNodes, msaOut);\r
 131         else\r
 132                 MakeRootMSA(v, NewTree, NewProgNodes, msaOut);\r
 133 \r
 134 #if     DEBUG\r
 135         AssertMSAEqIgnoreCaseAndGaps(msaIn, msaOut);\r
 136 #endif\r
 137 \r
 138         for (unsigned uNodeIndex = 0; uNodeIndex < uNodeCount; ++uNodeIndex)\r
 139                 DeleteProgNode(NewProgNodes[uNodeIndex]);\r
 140 \r
 141         delete[] NewProgNodes;\r
 142         }\r