+++ /dev/null
-#include "muscle.h"\r
-#include "msa.h"\r
-#include "tree.h"\r
-#include "clust.h"\r
-#include "profile.h"\r
-#include "pwpath.h"\r
-\r
-#define TRACE 0\r
-\r
-static void ProgressiveAlignSubfams(const Tree &tree, const unsigned Subfams[],\r
- unsigned uSubfamCount, const MSA SubfamMSAs[], MSA &msa);\r
-\r
-// Identify subfamilies in a tree.\r
-// Returns array of internal node indexes, one for each subfamily.\r
-// First try is to select groups by height (which should approximate\r
-// minimum percent identity), if this gives too many subfamilies then\r
-// we cut at a point that gives the maximum allowed number of subfams.\r
-static void GetSubfams(const Tree &tree, double dMaxHeight,\r
- unsigned uMaxSubfamCount, unsigned **ptrptrSubfams, unsigned *ptruSubfamCount)\r
- {\r
- const unsigned uNodeCount = tree.GetNodeCount();\r
-\r
- unsigned *Subfams = new unsigned[uNodeCount];\r
-\r
- unsigned uSubfamCount;\r
- ClusterByHeight(tree, dMaxHeight, Subfams, &uSubfamCount);\r
-\r
- if (uSubfamCount > uMaxSubfamCount)\r
- ClusterBySubfamCount(tree, uMaxSubfamCount, Subfams, &uSubfamCount);\r
-\r
- *ptrptrSubfams = Subfams;\r
- *ptruSubfamCount = uSubfamCount;\r
- }\r
-\r
-static void LogSubfams(const Tree &tree, const unsigned Subfams[],\r
- unsigned uSubfamCount)\r
- {\r
- const unsigned uNodeCount = tree.GetNodeCount();\r
- Log("%u subfamilies found\n", uSubfamCount);\r
- Log("Subfam Sequence\n");\r
- Log("------ --------\n");\r
- unsigned *Leaves = new unsigned[uNodeCount];\r
- for (unsigned uSubfamIndex = 0; uSubfamIndex < uSubfamCount; ++uSubfamIndex)\r
- {\r
- unsigned uSubfamNodeIndex = Subfams[uSubfamIndex];\r
- unsigned uLeafCount;\r
- GetLeaves(tree, uSubfamNodeIndex, Leaves, &uLeafCount);\r
- for (unsigned uLeafIndex = 0; uLeafIndex < uLeafCount; ++uLeafIndex)\r
- Log("%6u %s\n", uSubfamIndex + 1, tree.GetLeafName(Leaves[uLeafIndex]));\r
- Log("\n");\r
- }\r
- delete[] Leaves;\r
- }\r
-\r
-bool RefineSubfams(MSA &msa, const Tree &tree, unsigned uIters)\r
- {\r
- const unsigned uSeqCount = msa.GetSeqCount();\r
- if (uSeqCount < 3)\r
- return false;\r
-\r
- const double dMaxHeight = 0.6;\r
- const unsigned uMaxSubfamCount = 16;\r
- const unsigned uNodeCount = tree.GetNodeCount();\r
-\r
- unsigned *Subfams;\r
- unsigned uSubfamCount;\r
- GetSubfams(tree, dMaxHeight, uMaxSubfamCount, &Subfams, &uSubfamCount);\r
- assert(uSubfamCount <= uSeqCount);\r
-\r
- if (g_bVerbose)\r
- LogSubfams(tree, Subfams, uSubfamCount);\r
-\r
- MSA *SubfamMSAs = new MSA[uSubfamCount];\r
- unsigned *Leaves = new unsigned[uSeqCount];\r
- unsigned *Ids = new unsigned[uSeqCount];\r
-\r
- bool bAnyChanges = false;\r
- for (unsigned uSubfamIndex = 0; uSubfamIndex < uSubfamCount; ++uSubfamIndex)\r
- {\r
- unsigned uSubfam = Subfams[uSubfamIndex];\r
- unsigned uLeafCount;\r
- GetLeaves(tree, uSubfam, Leaves, &uLeafCount);\r
- assert(uLeafCount <= uSeqCount);\r
-\r
- LeafIndexesToIds(tree, Leaves, uLeafCount, Ids);\r
-\r
- MSA &msaSubfam = SubfamMSAs[uSubfamIndex];\r
- MSASubsetByIds(msa, Ids, uLeafCount, msaSubfam);\r
- DeleteGappedCols(msaSubfam);\r
-\r
-#if TRACE\r
- Log("Subfam %u MSA=\n", uSubfamIndex);\r
- msaSubfam.LogMe();\r
-#endif\r
-\r
- if (msaSubfam.GetSeqCount() <= 2)\r
- continue;\r
-\r
- // TODO /////////////////////////////////////////\r
- // Try using existing tree, may actually hurt to\r
- // re-estimate, may also be a waste of CPU & mem.\r
- /////////////////////////////////////////////////\r
- Tree SubfamTree;\r
- TreeFromMSA(msaSubfam, SubfamTree, g_Cluster2, g_Distance2, g_Root2);\r
-\r
- bool bAnyChangesThisSubfam;\r
- if (g_bAnchors)\r
- bAnyChangesThisSubfam = RefineVert(msaSubfam, SubfamTree, uIters);\r
- else\r
- bAnyChangesThisSubfam = RefineHoriz(msaSubfam, SubfamTree, uIters, false, false);\r
-#if TRACE\r
- Log("Subfam %u Changed %d\n", uSubfamIndex, bAnyChangesThisSubfam);\r
-#endif\r
- if (bAnyChangesThisSubfam)\r
- bAnyChanges = true;\r
- }\r
-\r
- if (bAnyChanges)\r
- ProgressiveAlignSubfams(tree, Subfams, uSubfamCount, SubfamMSAs, msa);\r
-\r
- delete[] Leaves;\r
- delete[] Subfams;\r
- delete[] SubfamMSAs;\r
-\r
- return bAnyChanges;\r
- }\r
-\r
-static void ProgressiveAlignSubfams(const Tree &tree, const unsigned Subfams[],\r
- unsigned uSubfamCount, const MSA SubfamMSAs[], MSA &msa)\r
- {\r
- const unsigned uNodeCount = tree.GetNodeCount();\r
-\r
- bool *Ready = new bool[uNodeCount];\r
- MSA **MSAs = new MSA *[uNodeCount];\r
- for (unsigned uNodeIndex = 0; uNodeIndex < uNodeCount; ++uNodeIndex)\r
- {\r
- Ready[uNodeIndex] = false;\r
- MSAs[uNodeIndex] = 0;\r
- }\r
-\r
- for (unsigned uSubfamIndex = 0; uSubfamIndex < uSubfamCount; ++uSubfamIndex)\r
- {\r
- unsigned uNodeIndex = Subfams[uSubfamIndex];\r
- Ready[uNodeIndex] = true;\r
- MSA *ptrMSA = new MSA;\r
- // TODO: Wasteful copy, needs re-design\r
- ptrMSA->Copy(SubfamMSAs[uSubfamIndex]);\r
- MSAs[uNodeIndex] = ptrMSA;\r
- }\r
-\r
- for (unsigned uNodeIndex = tree.FirstDepthFirstNode();\r
- NULL_NEIGHBOR != uNodeIndex;\r
- uNodeIndex = tree.NextDepthFirstNode(uNodeIndex))\r
- {\r
- if (tree.IsLeaf(uNodeIndex))\r
- continue;\r
-\r
- unsigned uRight = tree.GetRight(uNodeIndex);\r
- unsigned uLeft = tree.GetLeft(uNodeIndex);\r
- if (!Ready[uRight] || !Ready[uLeft])\r
- continue;\r
-\r
- MSA *ptrLeft = MSAs[uLeft];\r
- MSA *ptrRight = MSAs[uRight];\r
- assert(ptrLeft != 0 && ptrRight != 0);\r
-\r
- MSA *ptrParent = new MSA;\r
-\r
- PWPath Path;\r
- AlignTwoMSAs(*ptrLeft, *ptrRight, *ptrParent, Path);\r
-\r
- MSAs[uNodeIndex] = ptrParent;\r
- Ready[uNodeIndex] = true;\r
- Ready[uLeft] = false;\r
- Ready[uRight] = false;\r
-\r
- delete MSAs[uLeft];\r
- delete MSAs[uRight];\r
- MSAs[uLeft] = 0;\r
- MSAs[uRight] = 0;\r
- }\r
-\r
-#if DEBUG\r
- {\r
- unsigned uReadyCount = 0;\r
- for (unsigned uNodeIndex = 0; uNodeIndex < uNodeCount; ++uNodeIndex)\r
- {\r
- if (Ready[uNodeIndex])\r
- {\r
- assert(tree.IsRoot(uNodeIndex));\r
- ++uReadyCount;\r
- assert(0 != MSAs[uNodeIndex]);\r
- }\r
- else\r
- assert(0 == MSAs[uNodeIndex]);\r
- }\r
- assert(1 == uReadyCount);\r
- }\r
-#endif\r
-\r
- const unsigned uRoot = tree.GetRootNodeIndex();\r
- MSA *ptrRootAlignment = MSAs[uRoot];\r
-\r
- msa.Copy(*ptrRootAlignment);\r
-\r
- delete ptrRootAlignment;\r
-\r
-#if TRACE\r
- Log("After refine subfamilies, root alignment=\n");\r
- msa.LogMe();\r
-#endif\r
- }\r
git://source.jalview.org / jabaws.git / blobdiff