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