refinehoriz.cpp

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

#include "muscle.h"
#include "tree.h"
#include "msa.h"
#include "pwpath.h"
#include "profile.h"
#include "scorehistory.h"
#include "objscore.h"

unsigned g_uRefineHeightSubtree;
unsigned g_uRefineHeightSubtreeTotal;

#define TRACE			0
#define DIFFOBJSCORE	0

static bool TryRealign(MSA &msaIn, const Tree &tree, const unsigned Leaves1[],
  unsigned uCount1, const unsigned Leaves2[], unsigned uCount2,
  SCORE *ptrscoreBefore, SCORE *ptrscoreAfter,
  bool bLockLeft, bool bLockRight)
	{
#if	TRACE
	Log("TryRealign, msaIn=\n");
	msaIn.LogMe();
#endif

	const unsigned uSeqCount = msaIn.GetSeqCount();

	unsigned *Ids1 = new unsigned[uSeqCount];
	unsigned *Ids2 = new unsigned[uSeqCount];

	LeafIndexesToIds(tree, Leaves1, uCount1, Ids1);
	LeafIndexesToIds(tree, Leaves2, uCount2, Ids2);

	MSA msa1;
	MSA msa2;

	MSASubsetByIds(msaIn, Ids1, uCount1, msa1);
	MSASubsetByIds(msaIn, Ids2, uCount2, msa2);

#if	DEBUG
	ValidateMuscleIds(msa1);
	ValidateMuscleIds(msa2);
#endif

// Computing the objective score may be expensive for
// large numbers of sequences. As a speed optimization,
// we check whether the alignment changes. If it does
// not change, there is no need to compute the objective
// score. We test for the alignment changing by comparing
// the Viterbi paths before and after re-aligning.
	PWPath pathBefore;
	pathBefore.FromMSAPair(msa1, msa2);

	DeleteGappedCols(msa1);
	DeleteGappedCols(msa2);

	if (0 == msa1.GetColCount() || 0 == msa2.GetColCount())
		return false;

	MSA msaRealigned;
	PWPath pathAfter;

	AlignTwoMSAs(msa1, msa2, msaRealigned, pathAfter, bLockLeft, bLockRight);

	bool bAnyChanges = !pathAfter.Equal(pathBefore);
	unsigned uDiffCount1;
	unsigned uDiffCount2;
	static unsigned Edges1[10000];
	static unsigned Edges2[10000];
	DiffPaths(pathBefore, pathAfter, Edges1, &uDiffCount1, Edges2, &uDiffCount2);

#if	TRACE
	Log("TryRealign, msa1=\n");
	msa1.LogMe();
	Log("\nmsa2=\n");
	msa2.LogMe();
	Log("\nRealigned (changes %s)=\n", bAnyChanges ? "TRUE" : "FALSE");
	msaRealigned.LogMe();
#endif

	if (!bAnyChanges)
		{
		*ptrscoreBefore = 0;
		*ptrscoreAfter = 0;
		return false;
		}

	SetMSAWeightsMuscle(msaIn);
	SetMSAWeightsMuscle(msaRealigned);

#if	DIFFOBJSCORE
	const SCORE scoreDiff = DiffObjScore(msaIn, pathBefore, Edges1, uDiffCount1,
	  msaRealigned, pathAfter, Edges2, uDiffCount2);
	bool bAccept = (scoreDiff > 0);
	*ptrscoreBefore = 0;
	*ptrscoreAfter = scoreDiff;
	//const SCORE scoreBefore = ObjScoreIds(msaIn, Ids1, uCount1, Ids2, uCount2);
	//const SCORE scoreAfter = ObjScoreIds(msaRealigned, Ids1, uCount1, Ids2, uCount2);
	//Log("Diff = %.3g %.3g\n", scoreDiff, scoreAfter - scoreBefore);
#else
	const SCORE scoreBefore = ObjScoreIds(msaIn, Ids1, uCount1, Ids2, uCount2);
	const SCORE scoreAfter = ObjScoreIds(msaRealigned, Ids1, uCount1, Ids2, uCount2);

	bool bAccept = (scoreAfter > scoreBefore);

#if	TRACE
	Log("Score %g -> %g Accept %s\n", scoreBefore, scoreAfter, bAccept ? "TRUE" : "FALSE");
#endif

	*ptrscoreBefore = scoreBefore;
	*ptrscoreAfter = scoreAfter;
#endif

	if (bAccept)
		msaIn.Copy(msaRealigned);
	delete[] Ids1;
	delete[] Ids2;
	return bAccept;
	}

static void RefineHeightParts(MSA &msaIn, const Tree &tree,
 const unsigned InternalNodeIndexes[], bool bReversed, bool bRight,
 unsigned uIter, 
 ScoreHistory &History,
 bool *ptrbAnyChanges, bool *ptrbOscillating, bool bLockLeft, bool bLockRight)
	{
	*ptrbOscillating = false;

	const unsigned uSeqCount = msaIn.GetSeqCount();
	const unsigned uInternalNodeCount = uSeqCount - 1;

	unsigned *Leaves1 = new unsigned[uSeqCount];
	unsigned *Leaves2 = new unsigned[uSeqCount];

	const unsigned uRootNodeIndex = tree.GetRootNodeIndex();
	bool bAnyAccepted = false;
	for (unsigned i = 0; i < uInternalNodeCount; ++i)
		{
		const unsigned uInternalNodeIndex = InternalNodeIndexes[i];
		unsigned uNeighborNodeIndex;
		if (tree.IsRoot(uInternalNodeIndex) && !bRight)
			continue;
		else if (bRight)
			uNeighborNodeIndex = tree.GetRight(uInternalNodeIndex);
		else
			uNeighborNodeIndex = tree.GetLeft(uInternalNodeIndex);

		g_uTreeSplitNode1 = uInternalNodeIndex;
		g_uTreeSplitNode2 = uNeighborNodeIndex;

		unsigned uCount1;
		unsigned uCount2;

		GetLeaves(tree, uNeighborNodeIndex, Leaves1, &uCount1);
		GetLeavesExcluding(tree, uRootNodeIndex, uNeighborNodeIndex,
		  Leaves2, &uCount2);

#if	TRACE
		Log("\nRefineHeightParts node %u\n", uInternalNodeIndex);
		Log("Group1=");
		for (unsigned n = 0; n < uCount1; ++n)
			Log(" %u(%s)", Leaves1[n], tree.GetName(Leaves1[n]));
		Log("\n");
		Log("Group2=");
		for (unsigned n = 0; n < uCount2; ++n)
			Log(" %u(%s)", Leaves2[n], tree.GetName(Leaves2[n]));
		Log("\n");
#endif

		SCORE scoreBefore;
		SCORE scoreAfter;
		bool bAccepted = TryRealign(msaIn, tree, Leaves1, uCount1, Leaves2, uCount2,
		  &scoreBefore, &scoreAfter, bLockLeft, bLockRight);
		SetCurrentAlignment(msaIn);

		++g_uRefineHeightSubtree;
		Progress(g_uRefineHeightSubtree, g_uRefineHeightSubtreeTotal);

#if	TRACE
		if (uIter > 0)
			Log("Before %g %g\n", scoreBefore,
			  History.GetScore(uIter - 1, uInternalNodeIndex, bReversed, bRight));
#endif
		SCORE scoreMax = scoreAfter > scoreBefore? scoreAfter : scoreBefore;
		bool bRepeated = History.SetScore(uIter, uInternalNodeIndex, bRight, scoreMax);
		if (bRepeated)
			{
			*ptrbOscillating = true;
			break;
			}

		if (bAccepted)
			bAnyAccepted = true;
		}

	delete[] Leaves1;
	delete[] Leaves2;

	*ptrbAnyChanges = bAnyAccepted;
	}

// Return true if any changes made
bool RefineHoriz(MSA &msaIn, const Tree &tree, unsigned uIters, bool bLockLeft,
  bool bLockRight)
	{
#if	TRACE
	tree.LogMe();
#endif

	if (!tree.IsRooted())
		Quit("RefineHeight: requires rooted tree");

	const unsigned uSeqCount = msaIn.GetSeqCount();
	if (uSeqCount < 3)
		return false;

	const unsigned uInternalNodeCount = uSeqCount - 1;
	unsigned *InternalNodeIndexes = new unsigned[uInternalNodeCount];
	unsigned *InternalNodeIndexesR = new unsigned[uInternalNodeCount];

	GetInternalNodesInHeightOrder(tree, InternalNodeIndexes);

	ScoreHistory History(uIters, 2*uSeqCount - 1);

	bool bAnyChangesAnyIter = false;
	for (unsigned n = 0; n < uInternalNodeCount; ++n)
		InternalNodeIndexesR[uInternalNodeCount - 1 - n] = InternalNodeIndexes[n];

	for (unsigned uIter = 0; uIter < uIters; ++uIter)
		{
		bool bAnyChangesThisIter = false;
		IncIter();
		SetProgressDesc("Refine biparts");
		g_uRefineHeightSubtree = 0;
		g_uRefineHeightSubtreeTotal = uInternalNodeCount*2 - 1;

		bool bReverse = (uIter%2 != 0);
		unsigned *Internals;
		if (bReverse)
			Internals = InternalNodeIndexesR;
		else
			Internals = InternalNodeIndexes;

		bool bOscillating;
		for (unsigned i = 0; i < 2; ++i)
			{
			bool bAnyChanges = false;
			bool bRight;
			switch (i)
				{
			case 0:
				bRight = true;
				break;
			case 1:
				bRight = false;
				break;
			default:
				Quit("RefineHeight default case");
				}
			RefineHeightParts(msaIn, tree, Internals, bReverse, bRight,
			  uIter, 
			  History, 
			  &bAnyChanges, &bOscillating, bLockLeft, bLockRight);
			if (bOscillating)
				{
				ProgressStepsDone();
				goto Osc;
				}
			if (bAnyChanges)
				{
				bAnyChangesThisIter = true;
				bAnyChangesAnyIter = true;
				}
			}

		ProgressStepsDone();
		if (bOscillating)
			break;

		if (!bAnyChangesThisIter)
			break;
		}

Osc:
	delete[] InternalNodeIndexes;
	delete[] InternalNodeIndexesR;

	return bAnyChangesAnyIter;
	}