+++ /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