realigndiffse.cpp

unix,linux下编译。用于蛋白质

#include "muscle.h"
#include "msa.h"
#include "tree.h"
#include "profile.h"
#include "pwpath.h"
#include "seqvect.h"
#include "estring.h"

#define TRACE		0

void DeleteProgNode(ProgNode &Node)
	{
	delete[] Node.m_Prof;
	delete[] Node.m_EstringL;
	delete[] Node.m_EstringR;

	Node.m_Prof = 0;
	Node.m_EstringL = 0;
	Node.m_EstringR = 0;
	}

static void MakeNode(ProgNode &OldNode, ProgNode &NewNode, bool bSwapLR)
	{
	if (bSwapLR)
		{
		NewNode.m_EstringL = OldNode.m_EstringR;
		NewNode.m_EstringR = OldNode.m_EstringL;
		}
	else
		{
		NewNode.m_EstringL = OldNode.m_EstringL;
		NewNode.m_EstringR = OldNode.m_EstringR;
		}
	NewNode.m_Prof = OldNode.m_Prof;
	NewNode.m_uLength = OldNode.m_uLength;
	NewNode.m_Weight = OldNode.m_Weight;

	OldNode.m_Prof = 0;
	OldNode.m_EstringL = 0;
	OldNode.m_EstringR = 0;
	}

void RealignDiffsE(const MSA &msaIn, const SeqVect &v,
  const Tree &NewTree, const Tree &OldTree, 
  const unsigned uNewNodeIndexToOldNodeIndex[],
  MSA &msaOut, ProgNode *OldProgNodes)
	{
	assert(OldProgNodes != 0);

	const unsigned uNodeCount = NewTree.GetNodeCount();
	if (uNodeCount%2 == 0)
		Quit("RealignDiffs: Expected odd number of nodes");

	const unsigned uMergeCount = (uNodeCount - 1)/2;
	ProgNode *NewProgNodes = new ProgNode[uNodeCount];

	for (unsigned uNewNodeIndex = 0; uNewNodeIndex < uNodeCount; ++uNewNodeIndex)
		{
		if (NODE_CHANGED == uNewNodeIndexToOldNodeIndex[uNewNodeIndex])
			continue;

		unsigned uOldNodeIndex = uNewNodeIndexToOldNodeIndex[uNewNodeIndex];
		assert(uNewNodeIndex < uNodeCount);
		assert(uOldNodeIndex < uNodeCount);

		ProgNode &NewNode = NewProgNodes[uNewNodeIndex];
		ProgNode &OldNode = OldProgNodes[uOldNodeIndex];
		bool bSwapLR = false;
		if (!NewTree.IsLeaf(uNewNodeIndex))
			{
			unsigned uNewLeft = NewTree.GetLeft(uNewNodeIndex);
			unsigned uNewRight = NewTree.GetRight(uNewNodeIndex);
			unsigned uOld = uNewNodeIndexToOldNodeIndex[uNewNodeIndex];
			unsigned uOldLeft = OldTree.GetLeft(uOld);
			unsigned uOldRight = OldTree.GetRight(uOld);
			assert(uOldLeft < uNodeCount && uOldRight < uNodeCount);
			if (uOldLeft != uNewNodeIndexToOldNodeIndex[uNewLeft])
				{
				assert(uOldLeft == uNewNodeIndexToOldNodeIndex[uNewRight]);
				bSwapLR = true;
				}
			}
		MakeNode(OldNode, NewNode, bSwapLR);
#if	TRACE
		Log("MakeNode old=%u new=%u swap=%d length=%u weight=%.3g\n",
		  uOldNodeIndex, uNewNodeIndex, bSwapLR, NewNode.m_uLength, NewNode.m_Weight);
#endif
		}

	unsigned uJoin = 0;
	SetProgressDesc("Refine tree");
	for (unsigned uNewNodeIndex = NewTree.FirstDepthFirstNode();
	  NULL_NEIGHBOR != uNewNodeIndex;
	  uNewNodeIndex = NewTree.NextDepthFirstNode(uNewNodeIndex))
		{
		if (NODE_CHANGED != uNewNodeIndexToOldNodeIndex[uNewNodeIndex])
			continue;

		Progress(uJoin, uMergeCount - 1);
		++uJoin;

		const unsigned uMergeNodeIndex = uNewNodeIndex;
		ProgNode &Parent = NewProgNodes[uMergeNodeIndex];

		const unsigned uLeft = NewTree.GetLeft(uNewNodeIndex);
		const unsigned uRight = NewTree.GetRight(uNewNodeIndex);

		ProgNode &Node1 = NewProgNodes[uLeft];
		ProgNode &Node2 = NewProgNodes[uRight];

		AlignTwoProfs(
			Node1.m_Prof, Node1.m_uLength, Node1.m_Weight,
			Node2.m_Prof, Node2.m_uLength, Node2.m_Weight,
			Parent.m_Path,
			&Parent.m_Prof, &Parent.m_uLength);
		PathToEstrings(Parent.m_Path, &Parent.m_EstringL, &Parent.m_EstringR);

		Parent.m_Weight = Node1.m_Weight + Node2.m_Weight;

		delete[] Node1.m_Prof;
		delete[] Node2.m_Prof;

		Node1.m_Prof = 0;
		Node2.m_Prof = 0;
		}

	ProgressStepsDone();

	if (g_bBrenner)
		MakeRootMSABrenner((SeqVect &) v, NewTree, NewProgNodes, msaOut);
	else
		MakeRootMSA(v, NewTree, NewProgNodes, msaOut);

#if	DEBUG
	AssertMSAEqIgnoreCaseAndGaps(msaIn, msaOut);
#endif

	for (unsigned uNodeIndex = 0; uNodeIndex < uNodeCount; ++uNodeIndex)
		DeleteProgNode(NewProgNodes[uNodeIndex]);

	delete[] NewProgNodes;
	}