binaries/src/muscle/realigndiffs.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 \r
   7 #define TRACE   0\r
   8 \r
   9 // Progressive alignment according to a diffs tree.\r
  10 \r
  11 static void MakeNode(const MSA &msaIn, const Tree &Diffs, unsigned uDiffsNodeIndex,\r
  12    const unsigned IdToDiffsTreeNodeIndex[], ProgNode &Node)\r
  13         {\r
  14         const unsigned uSeqCount = msaIn.GetSeqCount();\r
  15 \r
  16         unsigned *Ids = new unsigned[uSeqCount];\r
  17 \r
  18         unsigned uSeqsInDiffCount = 0;\r
  19         for (unsigned uId = 0; uId < uSeqCount; ++uId)\r
  20                 {\r
  21                 if (IdToDiffsTreeNodeIndex[uId] == uDiffsNodeIndex)\r
  22                         {\r
  23                         Ids[uSeqsInDiffCount] = uId;\r
  24                         ++uSeqsInDiffCount;\r
  25                         }\r
  26                 }\r
  27         if (0 == uSeqsInDiffCount)\r
  28                 Quit("MakeNode: no seqs in diff");\r
  29 \r
  30         MSASubsetByIds(msaIn, Ids, uSeqsInDiffCount, Node.m_MSA);\r
  31 \r
  32 #if     DEBUG\r
  33         ValidateMuscleIds(Node.m_MSA);\r
  34 #endif\r
  35 \r
  36         DeleteGappedCols(Node.m_MSA);\r
  37         delete[] Ids;\r
  38         }\r
  39 \r
  40 void RealignDiffs(const MSA &msaIn, const Tree &Diffs,\r
  41   const unsigned IdToDiffsTreeNodeIndex[], MSA &msaOut)\r
  42         {\r
  43         assert(Diffs.IsRooted());\r
  44 \r
  45 #if     TRACE\r
  46         Log("RealignDiffs\n");\r
  47         Log("Diff tree:\n");\r
  48         Diffs.LogMe();\r
  49 #endif\r
  50 \r
  51         const unsigned uNodeCount = Diffs.GetNodeCount();\r
  52         if (uNodeCount%2 == 0)\r
  53                 Quit("RealignDiffs: Expected odd number of nodes");\r
  54 \r
  55         const unsigned uMergeCount = (uNodeCount - 1)/2;\r
  56 \r
  57         ProgNode *ProgNodes = new ProgNode[uNodeCount];\r
  58 \r
  59         unsigned uJoin = 0;\r
  60         SetProgressDesc("Refine tree");\r
  61         for (unsigned uDiffsNodeIndex = Diffs.FirstDepthFirstNode();\r
  62           NULL_NEIGHBOR != uDiffsNodeIndex;\r
  63           uDiffsNodeIndex = Diffs.NextDepthFirstNode(uDiffsNodeIndex))\r
  64                 {\r
  65                 if (Diffs.IsLeaf(uDiffsNodeIndex))\r
  66                         {\r
  67                         assert(uDiffsNodeIndex < uNodeCount);\r
  68                         if (uDiffsNodeIndex >= uNodeCount)\r
  69                                 Quit("TreeNodeIndex=%u NodeCount=%u\n", uDiffsNodeIndex, uNodeCount);\r
  70 \r
  71                         ProgNode &Node = ProgNodes[uDiffsNodeIndex];\r
  72                         MakeNode(msaIn, Diffs, uDiffsNodeIndex, IdToDiffsTreeNodeIndex, Node);\r
  73 \r
  74                         Node.m_uLength = Node.m_MSA.GetColCount();\r
  75                         }\r
  76                 else\r
  77                         {\r
  78                         Progress(uJoin, uMergeCount);\r
  79                         ++uJoin;\r
  80                         const unsigned uMergeNodeIndex = uDiffsNodeIndex;\r
  81                         ProgNode &Parent = ProgNodes[uMergeNodeIndex];\r
  82 \r
  83                         const unsigned uLeft = Diffs.GetLeft(uDiffsNodeIndex);\r
  84                         const unsigned uRight = Diffs.GetRight(uDiffsNodeIndex);\r
  85 \r
  86                         ProgNode &Node1 = ProgNodes[uLeft];\r
  87                         ProgNode &Node2 = ProgNodes[uRight];\r
  88 \r
  89                         PWPath Path;\r
  90                         AlignTwoMSAs(Node1.m_MSA, Node2.m_MSA, Parent.m_MSA, Path);\r
  91 \r
  92 #if     TRACE\r
  93                         {\r
  94                         Log("Combined:\n");\r
  95                         Parent.m_MSA.LogMe();\r
  96                         }\r
  97 #endif\r
  98 \r
  99                         Node1.m_MSA.Clear();\r
 100                         Node2.m_MSA.Clear();\r
 101                         }\r
 102                 }\r
 103         ProgressStepsDone();\r
 104 \r
 105         unsigned uRootNodeIndex = Diffs.GetRootNodeIndex();\r
 106         const ProgNode &RootProgNode = ProgNodes[uRootNodeIndex];\r
 107         msaOut.Copy(RootProgNode.m_MSA);\r
 108 \r
 109 #if     DEBUG\r
 110         AssertMSAEqIgnoreCaseAndGaps(msaIn, msaOut);\r
 111 #endif\r
 112 \r
 113         delete[] ProgNodes;\r
 114         ProgNodes = 0;\r
 115         }\r