--- /dev/null
+#include "muscle.h"\r
+#include "tree.h"\r
+#include "seqvect.h"\r
+#include "profile.h"\r
+#include "msa.h"\r
+#include "pwpath.h"\r
+#include "estring.h"\r
+\r
+#define TRACE 0\r
+#define VALIDATE 0\r
+\r
+static void PathSeq(const Seq &s, const PWPath &Path, bool bRight, Seq &sOut)\r
+ {\r
+ short *esA;\r
+ short *esB;\r
+ PathToEstrings(Path, &esA, &esB);\r
+\r
+ const unsigned uSeqLength = s.Length();\r
+ const unsigned uEdgeCount = Path.GetEdgeCount();\r
+\r
+ sOut.Clear();\r
+ sOut.SetName(s.GetName());\r
+ unsigned uPos = 0;\r
+ for (unsigned uEdgeIndex = 0; uEdgeIndex < uEdgeCount; ++uEdgeIndex)\r
+ {\r
+ const PWEdge &Edge = Path.GetEdge(uEdgeIndex);\r
+ char cType = Edge.cType;\r
+ if (bRight)\r
+ {\r
+ if (cType == 'I')\r
+ cType = 'D';\r
+ else if (cType == 'D')\r
+ cType = 'I';\r
+ }\r
+ switch (cType)\r
+ {\r
+ case 'M':\r
+ sOut.AppendChar(s[uPos++]);\r
+ break;\r
+ case 'D':\r
+ sOut.AppendChar('-');\r
+ break;\r
+ case 'I':\r
+ sOut.AppendChar(s[uPos++]);\r
+ break;\r
+ default:\r
+ Quit("PathSeq, invalid edge type %c", cType);\r
+ }\r
+ }\r
+ }\r
+\r
+#if VALIDATE\r
+\r
+static void MakeRootSeq(const Seq &s, const Tree &GuideTree, unsigned uLeafNodeIndex,\r
+ const ProgNode Nodes[], Seq &sRoot)\r
+ {\r
+ sRoot.Copy(s);\r
+ unsigned uNodeIndex = uLeafNodeIndex;\r
+ for (;;)\r
+ {\r
+ unsigned uParent = GuideTree.GetParent(uNodeIndex);\r
+ if (NULL_NEIGHBOR == uParent)\r
+ break;\r
+ bool bRight = (GuideTree.GetLeft(uParent) == uNodeIndex);\r
+ uNodeIndex = uParent;\r
+ const PWPath &Path = Nodes[uNodeIndex].m_Path;\r
+ Seq sTmp;\r
+ PathSeq(sRoot, Path, bRight, sTmp);\r
+ sTmp.SetId(0);\r
+ sRoot.Copy(sTmp);\r
+ }\r
+ }\r
+\r
+#endif // VALIDATE\r
+\r
+static short *MakeRootSeqE(const Seq &s, const Tree &GuideTree, unsigned uLeafNodeIndex,\r
+ const ProgNode Nodes[], Seq &sRoot, short *Estring1, short *Estring2)\r
+ {\r
+ short *EstringCurr = Estring1;\r
+ short *EstringNext = Estring2;\r
+\r
+ const unsigned uSeqLength = s.Length();\r
+ EstringCurr[0] = uSeqLength;\r
+ EstringCurr[1] = 0;\r
+\r
+ unsigned uNodeIndex = uLeafNodeIndex;\r
+ for (;;)\r
+ {\r
+ unsigned uParent = GuideTree.GetParent(uNodeIndex);\r
+ if (NULL_NEIGHBOR == uParent)\r
+ break;\r
+ bool bRight = (GuideTree.GetLeft(uParent) == uNodeIndex);\r
+ uNodeIndex = uParent;\r
+ const PWPath &Path = Nodes[uNodeIndex].m_Path;\r
+ const short *EstringNode = bRight ?\r
+ Nodes[uNodeIndex].m_EstringL : Nodes[uNodeIndex].m_EstringR;\r
+\r
+ MulEstrings(EstringCurr, EstringNode, EstringNext);\r
+#if TRACE\r
+ Log("\n");\r
+ Log("Curr=");\r
+ LogEstring(EstringCurr);\r
+ Log("\n");\r
+ Log("Node=");\r
+ LogEstring(EstringNode);\r
+ Log("\n");\r
+ Log("Prod=");\r
+ LogEstring(EstringNext);\r
+ Log("\n");\r
+#endif\r
+ short *EstringTmp = EstringNext;\r
+ EstringNext = EstringCurr;\r
+ EstringCurr = EstringTmp;\r
+ }\r
+ EstringOp(EstringCurr, s, sRoot);\r
+\r
+#if TRACE\r
+ Log("Root estring=");\r
+ LogEstring(EstringCurr);\r
+ Log("\n");\r
+ Log("Root seq=");\r
+ sRoot.LogMe();\r
+#endif\r
+ return EstringCurr;\r
+ }\r
+\r
+static unsigned GetFirstNodeIndex(const Tree &tree)\r
+ {\r
+ if (g_bStable)\r
+ return 0;\r
+ return tree.FirstDepthFirstNode();\r
+ }\r
+\r
+static unsigned GetNextNodeIndex(const Tree &tree, unsigned uPrevNodeIndex)\r
+ {\r
+ if (g_bStable)\r
+ {\r
+ const unsigned uNodeCount = tree.GetNodeCount();\r
+ unsigned uNodeIndex = uPrevNodeIndex;\r
+ for (;;)\r
+ {\r
+ ++uNodeIndex;\r
+ if (uNodeIndex >= uNodeCount)\r
+ return NULL_NEIGHBOR;\r
+ if (tree.IsLeaf(uNodeIndex))\r
+ return uNodeIndex;\r
+ }\r
+ }\r
+ unsigned uNodeIndex = uPrevNodeIndex;\r
+ for (;;)\r
+ {\r
+ uNodeIndex = tree.NextDepthFirstNode(uNodeIndex);\r
+ if (NULL_NEIGHBOR == uNodeIndex || tree.IsLeaf(uNodeIndex))\r
+ return uNodeIndex;\r
+ }\r
+ }\r
+\r
+void MakeRootMSA(const SeqVect &v, const Tree &GuideTree, ProgNode Nodes[],\r
+ MSA &a)\r
+ {\r
+#if TRACE\r
+ Log("MakeRootMSA Tree=");\r
+ GuideTree.LogMe();\r
+#endif\r
+ const unsigned uSeqCount = v.GetSeqCount();\r
+ unsigned uColCount = uInsane;\r
+ unsigned uSeqIndex = 0;\r
+ const unsigned uTreeNodeCount = GuideTree.GetNodeCount();\r
+ const unsigned uRootNodeIndex = GuideTree.GetRootNodeIndex();\r
+ const PWPath &RootPath = Nodes[uRootNodeIndex].m_Path;\r
+ const unsigned uRootColCount = RootPath.GetEdgeCount();\r
+ const unsigned uEstringSize = uRootColCount + 1;\r
+ short *Estring1 = new short[uEstringSize];\r
+ short *Estring2 = new short[uEstringSize];\r
+ SetProgressDesc("Root alignment");\r
+\r
+ unsigned uTreeNodeIndex = GetFirstNodeIndex(GuideTree);\r
+ do\r
+ {\r
+ Progress(uSeqIndex, uSeqCount);\r
+\r
+ unsigned uId = GuideTree.GetLeafId(uTreeNodeIndex);\r
+ const Seq &s = *(v[uId]);\r
+\r
+ Seq sRootE;\r
+ short *es = MakeRootSeqE(s, GuideTree, uTreeNodeIndex, Nodes, sRootE,\r
+ Estring1, Estring2);\r
+ Nodes[uTreeNodeIndex].m_EstringL = EstringNewCopy(es);\r
+\r
+#if VALIDATE\r
+ Seq sRoot;\r
+ MakeRootSeq(s, GuideTree, uTreeNodeIndex, Nodes, sRoot);\r
+ if (!sRoot.Eq(sRootE))\r
+ {\r
+ Log("sRoot=");\r
+ sRoot.LogMe();\r
+ Log("sRootE=");\r
+ sRootE.LogMe();\r
+ Quit("Root seqs differ");\r
+ }\r
+#if TRACE\r
+ Log("MakeRootSeq=\n");\r
+ sRoot.LogMe();\r
+#endif\r
+#endif\r
+\r
+ if (uInsane == uColCount)\r
+ {\r
+ uColCount = sRootE.Length();\r
+ a.SetSize(uSeqCount, uColCount);\r
+ }\r
+ else\r
+ {\r
+ assert(uColCount == sRootE.Length());\r
+ }\r
+ a.SetSeqName(uSeqIndex, s.GetName());\r
+ a.SetSeqId(uSeqIndex, uId);\r
+ for (unsigned uColIndex = 0; uColIndex < uColCount; ++uColIndex)\r
+ a.SetChar(uSeqIndex, uColIndex, sRootE[uColIndex]);\r
+ ++uSeqIndex;\r
+\r
+ uTreeNodeIndex = GetNextNodeIndex(GuideTree, uTreeNodeIndex);\r
+ }\r
+ while (NULL_NEIGHBOR != uTreeNodeIndex);\r
+\r
+ delete[] Estring1;\r
+ delete[] Estring2;\r
+\r
+ ProgressStepsDone();\r
+ assert(uSeqIndex == uSeqCount);\r
+ }\r
git://source.jalview.org / jabaws.git / blobdiff