sptest.cpp

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

#include "muscle.h"
#include "objscore.h"
#include "msa.h"
#include "textfile.h"
#include "pwpath.h"

const unsigned INDELS = 1;

static void GetPos(const char Str[], unsigned L, int *pi1, int *pi2)
	{
	int i1;
	for (;;)
		{
		i1 = rand()%(L-2) + 1;
		if (Str[i1] == 'M')
			break;
		}
	int i2;
	for (;;)
		{
		i2 = rand()%(L-2) + 1;
		if (i1 != i2 && Str[i2] == 'M')
			break;
		}
	*pi1 = i1;
	*pi2 = i2;
	}

static void MakePath(unsigned uSeqLength, unsigned uIndelCount, char Str[])
	{
	unsigned uPathLength = uSeqLength + uIndelCount;
	for (unsigned i = 0; i < uPathLength; ++i)
		Str[i] = 'M';

	for (unsigned i = 0; i < uIndelCount; ++i)
		{
		int i1, i2;
		GetPos(Str, uPathLength, &i1, &i2);
		Str[i1] = 'D';
		Str[i2] = 'I';
		}

	Str[uPathLength] = 0;
	Log("MakePath=%s\n", Str);
	}

void SPTest()
	{
	SetPPScore(PPSCORE_SV);

	SetListFileName("c:\\tmp\\muscle.log", false);

	TextFile file1("c:\\tmp\\msa1.afa");
	TextFile file2("c:\\tmp\\msa2.afa");

	MSA msa1;
	MSA msa2;

	msa1.FromFile(file1);
	msa2.FromFile(file2);

	Log("msa1=\n");
	msa1.LogMe();
	Log("msa2=\n");
	msa2.LogMe();

	const unsigned uColCount = msa1.GetColCount();
	if (msa2.GetColCount() != uColCount)
		Quit("Different lengths");

	const unsigned uSeqCount1 = msa1.GetSeqCount();
	const unsigned uSeqCount2 = msa2.GetSeqCount();
	const unsigned uSeqCount = uSeqCount1 + uSeqCount2;

	MSA::SetIdCount(uSeqCount);

	for (unsigned uSeqIndex1 = 0; uSeqIndex1 < uSeqCount1; ++uSeqIndex1)
		{
		msa1.SetSeqWeight(uSeqIndex1, 1.0);
		msa1.SetSeqId(uSeqIndex1, uSeqIndex1);
		}

	for (unsigned uSeqIndex2 = 0; uSeqIndex2 < uSeqCount2; ++uSeqIndex2)
		{
		msa2.SetSeqWeight(uSeqIndex2, 1.0);
		msa2.SetSeqId(uSeqIndex2, uSeqCount1 + uSeqIndex2);
		}

	MSA alnA;
	MSA alnB;

	char strPathA[1024];
	char strPathB[1024];
	MakePath(uColCount, INDELS, strPathA);
	MakePath(uColCount, INDELS, strPathB);

	PWPath PathA;
	PWPath PathB;
	PathA.FromStr(strPathA);
	PathB.FromStr(strPathB);

	Log("PathA=\n");
	PathA.LogMe();
	Log("PathB=\n");
	PathB.LogMe();

	AlignTwoMSAsGivenPath(PathA, msa1, msa2, alnA);
	AlignTwoMSAsGivenPath(PathB, msa1, msa2, alnB);

	for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)
		{
		alnA.SetSeqWeight(uSeqIndex, 1.0);
		alnB.SetSeqWeight(uSeqIndex, 1.0);
		}

	unsigned Seqs1[1024];
	unsigned Seqs2[1024];

	for (unsigned uSeqIndex1 = 0; uSeqIndex1 < uSeqCount1; ++uSeqIndex1)
		Seqs1[uSeqIndex1] = uSeqIndex1;

	for (unsigned uSeqIndex2 = 0; uSeqIndex2 < uSeqCount2; ++uSeqIndex2)
		Seqs2[uSeqIndex2] = uSeqCount1 + uSeqIndex2;

	MSA msaA1;
	MSA msaA2;
	MSA msaB1;
	MSA msaB2;
	MSAFromSeqSubset(alnA, Seqs1, uSeqCount1, msaA1);
	MSAFromSeqSubset(alnB, Seqs1, uSeqCount1, msaB1);
	MSAFromSeqSubset(alnA, Seqs2, uSeqCount2, msaA2);
	MSAFromSeqSubset(alnB, Seqs2, uSeqCount2, msaB2);

	for (unsigned uSeqIndex1 = 0; uSeqIndex1 < uSeqCount1; ++uSeqIndex1)
		{
		msaA1.SetSeqWeight(uSeqIndex1, 1.0);
		msaB1.SetSeqWeight(uSeqIndex1, 1.0);
		}

	for (unsigned uSeqIndex2 = 0; uSeqIndex2 < uSeqCount2; ++uSeqIndex2)
		{
		msaA2.SetSeqWeight(uSeqIndex2, 1.0);
		msaB2.SetSeqWeight(uSeqIndex2, 1.0);
		}

	Log("msaA1=\n");
	msaA1.LogMe();

	Log("msaB1=\n");
	msaB1.LogMe();

	Log("msaA2=\n");
	msaA2.LogMe();

	Log("msaB2=\n");
	msaB2.LogMe();

	Log("alnA=\n");
	alnA.LogMe();

	Log("AlnB=\n");
	alnB.LogMe();

	Log("\nSPA\n---\n");
	SCORE SPA = ObjScoreSP(alnA);
	Log("\nSPB\n---\n");
	SCORE SPB = ObjScoreSP(alnB);

	Log("\nXPA\n---\n");
	SCORE XPA = ObjScoreXP(msaA1, msaA2);
	Log("\nXPB\n---\n");
	SCORE XPB = ObjScoreXP(msaB1, msaB2);

	Log("SPA=%.4g SPB=%.4g Diff=%.4g\n", SPA, SPB, SPA - SPB);
	Log("XPA=%.4g XPB=%.4g Diff=%.4g\n", XPA, XPB, XPA - XPB);
	}