--- /dev/null
+#include "muscle.h"\r
+#include "tree.h"\r
+#include "textfile.h" // for test only\r
+#include "msa.h"\r
+#include "seqvect.h"\r
+#include "profile.h"\r
+#ifndef _MSC_VER\r
+#include <unistd.h> // for unlink\r
+#endif\r
+\r
+#define TRACE 0\r
+\r
+/***\r
+Find subfamilies from tree by following criteria:\r
+\r
+(a) number of leaves <= max,\r
+(b) is monophyletic, i.e. most recent common ancestor is parent\r
+of no more than one subfamily.\r
+***/\r
+\r
+static unsigned SubFamRecurse(const Tree &tree, unsigned uNodeIndex, unsigned uMaxLeafCount,\r
+ unsigned SubFams[], unsigned &uSubFamCount)\r
+ {\r
+ if (tree.IsLeaf(uNodeIndex))\r
+ return 1;\r
+\r
+ unsigned uLeft = tree.GetLeft(uNodeIndex);\r
+ unsigned uRight = tree.GetRight(uNodeIndex);\r
+ unsigned uLeftCount = SubFamRecurse(tree, uLeft, uMaxLeafCount, SubFams, uSubFamCount);\r
+ unsigned uRightCount = SubFamRecurse(tree, uRight, uMaxLeafCount, SubFams, uSubFamCount);\r
+\r
+ unsigned uLeafCount = uLeftCount + uRightCount;\r
+ if (uLeftCount + uRightCount > uMaxLeafCount)\r
+ {\r
+ if (uLeftCount <= uMaxLeafCount)\r
+ SubFams[uSubFamCount++] = uLeft;\r
+ if (uRightCount <= uMaxLeafCount)\r
+ SubFams[uSubFamCount++] = uRight;\r
+ }\r
+ else if (tree.IsRoot(uNodeIndex))\r
+ {\r
+ if (uSubFamCount != 0)\r
+ Quit("Error in SubFamRecurse");\r
+ SubFams[uSubFamCount++] = uNodeIndex;\r
+ }\r
+\r
+ return uLeafCount;\r
+ }\r
+\r
+void SubFam(const Tree &tree, unsigned uMaxLeafCount, unsigned SubFams[], unsigned *ptruSubFamCount)\r
+ {\r
+ *ptruSubFamCount = 0;\r
+ SubFamRecurse(tree, tree.GetRootNodeIndex(), uMaxLeafCount, SubFams, *ptruSubFamCount);\r
+\r
+#if TRACE\r
+ {\r
+ Log("\n");\r
+ Log("Tree:\n");\r
+ tree.LogMe();\r
+ //void DrawTree(const Tree &tree);\r
+ //DrawTree(tree);\r
+ Log("\n");\r
+ Log("%d subfams:\n", *ptruSubFamCount);\r
+ for (unsigned i = 0; i < *ptruSubFamCount; ++i)\r
+ Log(" %d=%d", i, SubFams[i]);\r
+ Log("\n");\r
+ }\r
+#endif\r
+ }\r
+\r
+//unsigned SubFams[9999];\r
+//unsigned uSubFamCount;\r
+//\r
+//static unsigned DistFromRoot(const Tree &tree, unsigned uNodeIndex)\r
+// {\r
+// const unsigned uRoot = tree.GetRootNodeIndex();\r
+// unsigned uDist = 0;\r
+// while (uNodeIndex != uRoot)\r
+// {\r
+// ++uDist;\r
+// uNodeIndex = tree.GetParent(uNodeIndex);\r
+// }\r
+// return uDist;\r
+// }\r
+//\r
+//static void DrawNode(const Tree &tree, unsigned uNodeIndex)\r
+// {\r
+// if (!tree.IsLeaf(uNodeIndex))\r
+// DrawNode(tree, tree.GetLeft(uNodeIndex));\r
+//\r
+// unsigned uDist = DistFromRoot(tree, uNodeIndex);\r
+// for (unsigned i = 0; i < 5*uDist; ++i)\r
+// Log(" ");\r
+// Log("%d", uNodeIndex);\r
+// for (unsigned i = 0; i < uSubFamCount; ++i)\r
+// if (uNodeIndex == SubFams[i])\r
+// {\r
+// Log("*");\r
+// break;\r
+// }\r
+// Log("\n");\r
+//\r
+// if (!tree.IsLeaf(uNodeIndex))\r
+// DrawNode(tree, tree.GetRight(uNodeIndex));\r
+// }\r
+//\r
+//static void DrawTree(const Tree &tree)\r
+// {\r
+// unsigned uRoot = tree.GetRootNodeIndex();\r
+// DrawNode(tree, uRoot);\r
+// }\r
+//\r
+//void TestSubFams(const char *FileName)\r
+// {\r
+// Tree tree;\r
+// TextFile f(FileName);\r
+// tree.FromFile(f);\r
+// SubFam(tree, 5, SubFams, &uSubFamCount);\r
+// DrawTree(tree);\r
+// }\r
+\r
+static void SetInFam(const Tree &tree, unsigned uNodeIndex, bool NodeInSubFam[])\r
+ {\r
+ if (tree.IsLeaf(uNodeIndex))\r
+ return;\r
+ unsigned uLeft = tree.GetLeft(uNodeIndex);\r
+ unsigned uRight = tree.GetRight(uNodeIndex);\r
+ NodeInSubFam[uLeft] = true;\r
+ NodeInSubFam[uRight] = true;\r
+\r
+ SetInFam(tree, uLeft, NodeInSubFam);\r
+ SetInFam(tree, uRight, NodeInSubFam);\r
+ }\r
+\r
+void AlignSubFam(SeqVect &vAll, const Tree &GuideTree, unsigned uNodeIndex,\r
+ MSA &msaOut)\r
+ {\r
+ const unsigned uSeqCount = vAll.GetSeqCount();\r
+\r
+ const char *InTmp = "asf_in.tmp";\r
+ const char *OutTmp = "asf_out.tmp";\r
+\r
+ unsigned *Leaves = new unsigned[uSeqCount];\r
+ unsigned uLeafCount;\r
+ GetLeaves(GuideTree, uNodeIndex, Leaves, &uLeafCount);\r
+\r
+ SeqVect v;\r
+ for (unsigned i = 0; i < uLeafCount; ++i)\r
+ {\r
+ unsigned uLeafNodeIndex = Leaves[i];\r
+ unsigned uId = GuideTree.GetLeafId(uLeafNodeIndex);\r
+ Seq &s = vAll.GetSeqById(uId);\r
+ v.AppendSeq(s);\r
+ }\r
+\r
+#if TRACE\r
+ {\r
+ Log("Align subfam[node=%d, size=%d] ", uNodeIndex, uLeafCount);\r
+ for (unsigned i = 0; i < uLeafCount; ++i)\r
+ Log(" %s", v.GetSeqName(i));\r
+ Log("\n");\r
+ }\r
+#endif\r
+\r
+ TextFile fIn(InTmp, true);\r
+\r
+ v.ToFASTAFile(fIn);\r
+ fIn.Close();\r
+\r
+ char CmdLine[4096];\r
+ sprintf(CmdLine, "probcons %s > %s 2> /dev/null", InTmp, OutTmp);\r
+// sprintf(CmdLine, "muscle -in %s -out %s -maxiters 1", InTmp, OutTmp);\r
+ int NotUsed = system(CmdLine);\r
+\r
+ TextFile fOut(OutTmp);\r
+ msaOut.FromFile(fOut);\r
+\r
+ for (unsigned uSeqIndex = 0; uSeqIndex < uLeafCount; ++uSeqIndex)\r
+ {\r
+ const char *Name = msaOut.GetSeqName(uSeqIndex);\r
+ unsigned uId = vAll.GetSeqIdFromName(Name);\r
+ msaOut.SetSeqId(uSeqIndex, uId);\r
+ }\r
+\r
+ unlink(InTmp);\r
+ unlink(OutTmp);\r
+\r
+ delete[] Leaves;\r
+ }\r
+\r
+void ProgAlignSubFams()\r
+ {\r
+ MSA msaOut;\r
+\r
+ SetOutputFileName(g_pstrOutFileName);\r
+ SetInputFileName(g_pstrInFileName);\r
+\r
+ SetMaxIters(g_uMaxIters);\r
+ SetSeqWeightMethod(g_SeqWeight1);\r
+\r
+ TextFile fileIn(g_pstrInFileName);\r
+ SeqVect v;\r
+ v.FromFASTAFile(fileIn);\r
+ const unsigned uSeqCount = v.Length();\r
+\r
+ if (0 == uSeqCount)\r
+ Quit("No sequences in input file");\r
+\r
+ ALPHA Alpha = ALPHA_Undefined;\r
+ switch (g_SeqType)\r
+ {\r
+ case SEQTYPE_Auto:\r
+ Alpha = v.GuessAlpha();\r
+ break;\r
+\r
+ case SEQTYPE_Protein:\r
+ Alpha = ALPHA_Amino;\r
+ break;\r
+\r
+ case SEQTYPE_DNA:\r
+ Alpha = ALPHA_DNA;\r
+ break;\r
+\r
+ case SEQTYPE_RNA:\r
+ Alpha = ALPHA_RNA;\r
+ break;\r
+\r
+ default:\r
+ Quit("Invalid seq type");\r
+ }\r
+ SetAlpha(Alpha);\r
+ v.FixAlpha();\r
+\r
+ PTR_SCOREMATRIX UserMatrix = 0;\r
+ if (0 != g_pstrMatrixFileName)\r
+ {\r
+ const char *FileName = g_pstrMatrixFileName;\r
+ const char *Path = getenv("MUSCLE_MXPATH");\r
+ if (Path != 0)\r
+ {\r
+ size_t n = strlen(Path) + 1 + strlen(FileName) + 1;\r
+ char *NewFileName = new char[n];\r
+ sprintf(NewFileName, "%s/%s", Path, FileName);\r
+ FileName = NewFileName;\r
+ }\r
+ TextFile File(FileName);\r
+ UserMatrix = ReadMx(File);\r
+ g_Alpha = ALPHA_Amino;\r
+ g_PPScore = PPSCORE_SP;\r
+ }\r
+\r
+ SetPPScore();\r
+\r
+ if (0 != UserMatrix)\r
+ g_ptrScoreMatrix = UserMatrix;\r
+\r
+ if (ALPHA_DNA == Alpha || ALPHA_RNA == Alpha)\r
+ {\r
+ SetPPScore(PPSCORE_SPN);\r
+ g_Distance1 = DISTANCE_Kmer4_6;\r
+ }\r
+\r
+ unsigned uMaxL = 0;\r
+ unsigned uTotL = 0;\r
+ for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
+ {\r
+ unsigned L = v.GetSeq(uSeqIndex).Length();\r
+ uTotL += L;\r
+ if (L > uMaxL)\r
+ uMaxL = L;\r
+ }\r
+\r
+ SetIter(1);\r
+ g_bDiags = g_bDiags1;\r
+ SetSeqStats(uSeqCount, uMaxL, uTotL/uSeqCount);\r
+\r
+ SetMuscleSeqVect(v);\r
+\r
+ MSA::SetIdCount(uSeqCount);\r
+\r
+// Initialize sequence ids.\r
+// From this point on, ids must somehow propogate from here.\r
+ for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
+ v.SetSeqId(uSeqIndex, uSeqIndex);\r
+\r
+ if (uSeqCount > 1)\r
+ MHackStart(v);\r
+\r
+ if (0 == uSeqCount)\r
+ {\r
+ msaOut.Clear();\r
+ return;\r
+ }\r
+\r
+ if (1 == uSeqCount && ALPHA_Amino == Alpha)\r
+ {\r
+ const Seq &s = v.GetSeq(0);\r
+ msaOut.FromSeq(s);\r
+ return;\r
+ }\r
+\r
+ Tree GuideTree;\r
+ TreeFromSeqVect(v, GuideTree, g_Cluster1, g_Distance1, g_Root1);\r
+ SetMuscleTree(GuideTree);\r
+\r
+ MSA msa;\r
+ if (g_bLow)\r
+ {\r
+ ProgNode *ProgNodes = 0;\r
+ ProgNodes = ProgressiveAlignE(v, GuideTree, msa);\r
+ delete[] ProgNodes;\r
+ }\r
+ else\r
+ ProgressiveAlign(v, GuideTree, msa);\r
+ SetCurrentAlignment(msa);\r
+ TreeFromMSA(msa, GuideTree, g_Cluster2, g_Distance2, g_Root2);\r
+ SetMuscleTree(GuideTree);\r
+\r
+ unsigned *SubFams = new unsigned[uSeqCount];\r
+ unsigned uSubFamCount;\r
+ SubFam(GuideTree, g_uMaxSubFamCount, SubFams, &uSubFamCount);\r
+\r
+ SetProgressDesc("Align node");\r
+ const unsigned uNodeCount = 2*uSeqCount - 1;\r
+\r
+ ProgNode *ProgNodes = new ProgNode[uNodeCount];\r
+ bool *NodeIsSubFam = new bool[uNodeCount];\r
+ bool *NodeInSubFam = new bool[uNodeCount];\r
+\r
+ for (unsigned i = 0; i < uNodeCount; ++i)\r
+ {\r
+ NodeIsSubFam[i] = false;\r
+ NodeInSubFam[i] = false;\r
+ }\r
+\r
+ for (unsigned i = 0; i < uSubFamCount; ++i)\r
+ {\r
+ unsigned uNodeIndex = SubFams[i];\r
+ assert(uNodeIndex < uNodeCount);\r
+ NodeIsSubFam[uNodeIndex] = true;\r
+ SetInFam(GuideTree, uNodeIndex, NodeInSubFam);\r
+ }\r
+\r
+ unsigned uJoin = 0;\r
+ unsigned uTreeNodeIndex = GuideTree.FirstDepthFirstNode();\r
+ do\r
+ {\r
+ if (NodeIsSubFam[uTreeNodeIndex])\r
+ {\r
+#if TRACE\r
+ Log("Node %d: align subfam\n", uTreeNodeIndex);\r
+#endif\r
+ ProgNode &Node = ProgNodes[uTreeNodeIndex];\r
+ AlignSubFam(v, GuideTree, uTreeNodeIndex, Node.m_MSA);\r
+ Node.m_uLength = Node.m_MSA.GetColCount();\r
+ }\r
+ else if (!NodeInSubFam[uTreeNodeIndex])\r
+ {\r
+#if TRACE\r
+ Log("Node %d: align two subfams\n", uTreeNodeIndex);\r
+#endif\r
+ Progress(uJoin, uSubFamCount - 1);\r
+ ++uJoin;\r
+\r
+ const unsigned uMergeNodeIndex = uTreeNodeIndex;\r
+ ProgNode &Parent = ProgNodes[uMergeNodeIndex];\r
+\r
+ const unsigned uLeft = GuideTree.GetLeft(uTreeNodeIndex);\r
+ const unsigned uRight = GuideTree.GetRight(uTreeNodeIndex);\r
+\r
+ ProgNode &Node1 = ProgNodes[uLeft];\r
+ ProgNode &Node2 = ProgNodes[uRight];\r
+\r
+ PWPath Path;\r
+ AlignTwoMSAs(Node1.m_MSA, Node2.m_MSA, Parent.m_MSA, Path);\r
+ Parent.m_uLength = Parent.m_MSA.GetColCount();\r
+\r
+ Node1.m_MSA.Clear();\r
+ Node2.m_MSA.Clear();\r
+ }\r
+ else\r
+ {\r
+#if TRACE\r
+ Log("Node %d: in subfam\n", uTreeNodeIndex);\r
+#endif\r
+ ;\r
+ }\r
+ uTreeNodeIndex = GuideTree.NextDepthFirstNode(uTreeNodeIndex);\r
+ }\r
+ while (NULL_NEIGHBOR != uTreeNodeIndex);\r
+ ProgressStepsDone();\r
+\r
+ unsigned uRootNodeIndex = GuideTree.GetRootNodeIndex();\r
+ ProgNode &RootProgNode = ProgNodes[uRootNodeIndex];\r
+\r
+ TextFile fOut(g_pstrOutFileName, true);\r
+ MHackEnd(RootProgNode.m_MSA);\r
+ RootProgNode.m_MSA.ToFile(fOut);\r
+\r
+ delete[] NodeInSubFam;\r
+ delete[] NodeIsSubFam;\r
+ delete[] ProgNodes;\r
+ delete[] SubFams;\r
+\r
+ ProgNodes = 0;\r
+ NodeInSubFam = 0;\r
+ NodeIsSubFam = 0;\r
+ SubFams = 0;\r
+ }\r
git://source.jalview.org / jabaws.git / blobdiff