website/archive/binaries/mac/src/muscle/refinehoriz.cpp

   1 #include "muscle.h"\r
   2 #include "tree.h"\r
   3 #include "msa.h"\r
   4 #include "pwpath.h"\r
   5 #include "profile.h"\r
   6 #include "scorehistory.h"\r
   7 #include "objscore.h"\r
   8 \r
   9 unsigned g_uRefineHeightSubtree;\r
  10 unsigned g_uRefineHeightSubtreeTotal;\r
  11 \r
  12 #define TRACE                   0\r
  13 #define DIFFOBJSCORE    0\r
  14 \r
  15 static bool TryRealign(MSA &msaIn, const Tree &tree, const unsigned Leaves1[],\r
  16   unsigned uCount1, const unsigned Leaves2[], unsigned uCount2,\r
  17   SCORE *ptrscoreBefore, SCORE *ptrscoreAfter,\r
  18   bool bLockLeft, bool bLockRight)\r
  19         {\r
  20 #if     TRACE\r
  21         Log("TryRealign, msaIn=\n");\r
  22         msaIn.LogMe();\r
  23 #endif\r
  24 \r
  25         const unsigned uSeqCount = msaIn.GetSeqCount();\r
  26 \r
  27         unsigned *Ids1 = new unsigned[uSeqCount];\r
  28         unsigned *Ids2 = new unsigned[uSeqCount];\r
  29 \r
  30         LeafIndexesToIds(tree, Leaves1, uCount1, Ids1);\r
  31         LeafIndexesToIds(tree, Leaves2, uCount2, Ids2);\r
  32 \r
  33         MSA msa1;\r
  34         MSA msa2;\r
  35 \r
  36         MSASubsetByIds(msaIn, Ids1, uCount1, msa1);\r
  37         MSASubsetByIds(msaIn, Ids2, uCount2, msa2);\r
  38 \r
  39 #if     DEBUG\r
  40         ValidateMuscleIds(msa1);\r
  41         ValidateMuscleIds(msa2);\r
  42 #endif\r
  43 \r
  44 // Computing the objective score may be expensive for\r
  45 // large numbers of sequences. As a speed optimization,\r
  46 // we check whether the alignment changes. If it does\r
  47 // not change, there is no need to compute the objective\r
  48 // score. We test for the alignment changing by comparing\r
  49 // the Viterbi paths before and after re-aligning.\r
  50         PWPath pathBefore;\r
  51         pathBefore.FromMSAPair(msa1, msa2);\r
  52 \r
  53         DeleteGappedCols(msa1);\r
  54         DeleteGappedCols(msa2);\r
  55 \r
  56         if (0 == msa1.GetColCount() || 0 == msa2.GetColCount())\r
  57                 return false;\r
  58 \r
  59         MSA msaRealigned;\r
  60         PWPath pathAfter;\r
  61 \r
  62         AlignTwoMSAs(msa1, msa2, msaRealigned, pathAfter, bLockLeft, bLockRight);\r
  63 \r
  64         bool bAnyChanges = !pathAfter.Equal(pathBefore);\r
  65         unsigned uDiffCount1;\r
  66         unsigned uDiffCount2;\r
  67         static unsigned Edges1[10000];\r
  68         static unsigned Edges2[10000];\r
  69         DiffPaths(pathBefore, pathAfter, Edges1, &uDiffCount1, Edges2, &uDiffCount2);\r
  70 \r
  71 #if     TRACE\r
  72         Log("TryRealign, msa1=\n");\r
  73         msa1.LogMe();\r
  74         Log("\nmsa2=\n");\r
  75         msa2.LogMe();\r
  76         Log("\nRealigned (changes %s)=\n", bAnyChanges ? "TRUE" : "FALSE");\r
  77         msaRealigned.LogMe();\r
  78 #endif\r
  79 \r
  80         if (!bAnyChanges)\r
  81                 {\r
  82                 *ptrscoreBefore = 0;\r
  83                 *ptrscoreAfter = 0;\r
  84                 return false;\r
  85                 }\r
  86 \r
  87         SetMSAWeightsMuscle(msaIn);\r
  88         SetMSAWeightsMuscle(msaRealigned);\r
  89 \r
  90 #if     DIFFOBJSCORE\r
  91         const SCORE scoreDiff = DiffObjScore(msaIn, pathBefore, Edges1, uDiffCount1,\r
  92           msaRealigned, pathAfter, Edges2, uDiffCount2);\r
  93         bool bAccept = (scoreDiff > 0);\r
  94         *ptrscoreBefore = 0;\r
  95         *ptrscoreAfter = scoreDiff;\r
  96         //const SCORE scoreBefore = ObjScoreIds(msaIn, Ids1, uCount1, Ids2, uCount2);\r
  97         //const SCORE scoreAfter = ObjScoreIds(msaRealigned, Ids1, uCount1, Ids2, uCount2);\r
  98         //Log("Diff = %.3g %.3g\n", scoreDiff, scoreAfter - scoreBefore);\r
  99 #else\r
 100         const SCORE scoreBefore = ObjScoreIds(msaIn, Ids1, uCount1, Ids2, uCount2);\r
 101         const SCORE scoreAfter = ObjScoreIds(msaRealigned, Ids1, uCount1, Ids2, uCount2);\r
 102 \r
 103         bool bAccept = (scoreAfter > scoreBefore);\r
 104 \r
 105 #if     TRACE\r
 106         Log("Score %g -> %g Accept %s\n", scoreBefore, scoreAfter, bAccept ? "TRUE" : "FALSE");\r
 107 #endif\r
 108 \r
 109         *ptrscoreBefore = scoreBefore;\r
 110         *ptrscoreAfter = scoreAfter;\r
 111 #endif\r
 112 \r
 113         if (bAccept)\r
 114                 msaIn.Copy(msaRealigned);\r
 115         delete[] Ids1;\r
 116         delete[] Ids2;\r
 117         return bAccept;\r
 118         }\r
 119 \r
 120 static void RefineHeightParts(MSA &msaIn, const Tree &tree,\r
 121  const unsigned InternalNodeIndexes[], bool bReversed, bool bRight,\r
 122  unsigned uIter, \r
 123  ScoreHistory &History,\r
 124  bool *ptrbAnyChanges, bool *ptrbOscillating, bool bLockLeft, bool bLockRight)\r
 125         {\r
 126         *ptrbOscillating = false;\r
 127 \r
 128         const unsigned uSeqCount = msaIn.GetSeqCount();\r
 129         const unsigned uInternalNodeCount = uSeqCount - 1;\r
 130 \r
 131         unsigned *Leaves1 = new unsigned[uSeqCount];\r
 132         unsigned *Leaves2 = new unsigned[uSeqCount];\r
 133 \r
 134         const unsigned uRootNodeIndex = tree.GetRootNodeIndex();\r
 135         bool bAnyAccepted = false;\r
 136         for (unsigned i = 0; i < uInternalNodeCount; ++i)\r
 137                 {\r
 138                 const unsigned uInternalNodeIndex = InternalNodeIndexes[i];\r
 139                 unsigned uNeighborNodeIndex;\r
 140                 if (tree.IsRoot(uInternalNodeIndex) && !bRight)\r
 141                         continue;\r
 142                 else if (bRight)\r
 143                         uNeighborNodeIndex = tree.GetRight(uInternalNodeIndex);\r
 144                 else\r
 145                         uNeighborNodeIndex = tree.GetLeft(uInternalNodeIndex);\r
 146 \r
 147                 g_uTreeSplitNode1 = uInternalNodeIndex;\r
 148                 g_uTreeSplitNode2 = uNeighborNodeIndex;\r
 149 \r
 150                 unsigned uCount1;\r
 151                 unsigned uCount2;\r
 152 \r
 153                 GetLeaves(tree, uNeighborNodeIndex, Leaves1, &uCount1);\r
 154                 GetLeavesExcluding(tree, uRootNodeIndex, uNeighborNodeIndex,\r
 155                   Leaves2, &uCount2);\r
 156 \r
 157 #if     TRACE\r
 158                 Log("\nRefineHeightParts node %u\n", uInternalNodeIndex);\r
 159                 Log("Group1=");\r
 160                 for (unsigned n = 0; n < uCount1; ++n)\r
 161                         Log(" %u(%s)", Leaves1[n], tree.GetName(Leaves1[n]));\r
 162                 Log("\n");\r
 163                 Log("Group2=");\r
 164                 for (unsigned n = 0; n < uCount2; ++n)\r
 165                         Log(" %u(%s)", Leaves2[n], tree.GetName(Leaves2[n]));\r
 166                 Log("\n");\r
 167 #endif\r
 168 \r
 169                 SCORE scoreBefore;\r
 170                 SCORE scoreAfter;\r
 171                 bool bAccepted = TryRealign(msaIn, tree, Leaves1, uCount1, Leaves2, uCount2,\r
 172                   &scoreBefore, &scoreAfter, bLockLeft, bLockRight);\r
 173                 SetCurrentAlignment(msaIn);\r
 174 \r
 175                 ++g_uRefineHeightSubtree;\r
 176                 Progress(g_uRefineHeightSubtree, g_uRefineHeightSubtreeTotal);\r
 177 \r
 178 #if     TRACE\r
 179                 if (uIter > 0)\r
 180                         Log("Before %g %g\n", scoreBefore,\r
 181                           History.GetScore(uIter - 1, uInternalNodeIndex, bReversed, bRight));\r
 182 #endif\r
 183                 SCORE scoreMax = scoreAfter > scoreBefore? scoreAfter : scoreBefore;\r
 184                 bool bRepeated = History.SetScore(uIter, uInternalNodeIndex, bRight, scoreMax);\r
 185                 if (bRepeated)\r
 186                         {\r
 187                         *ptrbOscillating = true;\r
 188                         break;\r
 189                         }\r
 190 \r
 191                 if (bAccepted)\r
 192                         bAnyAccepted = true;\r
 193                 }\r
 194 \r
 195         delete[] Leaves1;\r
 196         delete[] Leaves2;\r
 197 \r
 198         *ptrbAnyChanges = bAnyAccepted;\r
 199         }\r
 200 \r
 201 // Return true if any changes made\r
 202 bool RefineHoriz(MSA &msaIn, const Tree &tree, unsigned uIters, bool bLockLeft,\r
 203   bool bLockRight)\r
 204         {\r
 205 #if     TRACE\r
 206         tree.LogMe();\r
 207 #endif\r
 208 \r
 209         if (!tree.IsRooted())\r
 210                 Quit("RefineHeight: requires rooted tree");\r
 211 \r
 212         const unsigned uSeqCount = msaIn.GetSeqCount();\r
 213         if (uSeqCount < 3)\r
 214                 return false;\r
 215 \r
 216         const unsigned uInternalNodeCount = uSeqCount - 1;\r
 217         unsigned *InternalNodeIndexes = new unsigned[uInternalNodeCount];\r
 218         unsigned *InternalNodeIndexesR = new unsigned[uInternalNodeCount];\r
 219 \r
 220         GetInternalNodesInHeightOrder(tree, InternalNodeIndexes);\r
 221 \r
 222         ScoreHistory History(uIters, 2*uSeqCount - 1);\r
 223 \r
 224         bool bAnyChangesAnyIter = false;\r
 225         for (unsigned n = 0; n < uInternalNodeCount; ++n)\r
 226                 InternalNodeIndexesR[uInternalNodeCount - 1 - n] = InternalNodeIndexes[n];\r
 227 \r
 228         for (unsigned uIter = 0; uIter < uIters; ++uIter)\r
 229                 {\r
 230                 bool bAnyChangesThisIter = false;\r
 231                 IncIter();\r
 232                 SetProgressDesc("Refine biparts");\r
 233                 g_uRefineHeightSubtree = 0;\r
 234                 g_uRefineHeightSubtreeTotal = uInternalNodeCount*2 - 1;\r
 235 \r
 236                 bool bReverse = (uIter%2 != 0);\r
 237                 unsigned *Internals;\r
 238                 if (bReverse)\r
 239                         Internals = InternalNodeIndexesR;\r
 240                 else\r
 241                         Internals = InternalNodeIndexes;\r
 242 \r
 243                 bool bOscillating;\r
 244                 for (unsigned i = 0; i < 2; ++i)\r
 245                         {\r
 246                         bool bAnyChanges = false;\r
 247                         bool bRight;\r
 248                         switch (i)\r
 249                                 {\r
 250                         case 0:\r
 251                                 bRight = true;\r
 252                                 break;\r
 253                         case 1:\r
 254                                 bRight = false;\r
 255                                 break;\r
 256                         default:\r
 257                                 Quit("RefineHeight default case");\r
 258                                 }\r
 259                         RefineHeightParts(msaIn, tree, Internals, bReverse, bRight,\r
 260                           uIter, \r
 261                           History, \r
 262                           &bAnyChanges, &bOscillating, bLockLeft, bLockRight);\r
 263                         if (bOscillating)\r
 264                                 {\r
 265                                 ProgressStepsDone();\r
 266                                 goto Osc;\r
 267                                 }\r
 268                         if (bAnyChanges)\r
 269                                 {\r
 270                                 bAnyChangesThisIter = true;\r
 271                                 bAnyChangesAnyIter = true;\r
 272                                 }\r
 273                         }\r
 274 \r
 275                 ProgressStepsDone();\r
 276                 if (bOscillating)\r
 277                         break;\r
 278 \r
 279                 if (!bAnyChangesThisIter)\r
 280                         break;\r
 281                 }\r
 282 \r
 283 Osc:\r
 284         delete[] InternalNodeIndexes;\r
 285         delete[] InternalNodeIndexesR;\r
 286 \r
 287         return bAnyChangesAnyIter;\r
 288         }\r