⭐ 欢迎来到虫虫下载站！ | 📦 资源下载 📁 资源专辑 ℹ️ 关于我们
⭐ 虫虫下载站
📤 上传资源

📄 alignmentclass.cpp

📁 利用动态规划算法进行核酸序列的两两比对的源代码。
💻 CPP
字号:
🔍 
// aligmentClass.cpp: implementation of the AligmentClass class.
//
//////////////////////////////////////////////////////////////////////


#include "math.h"
#include "AlignmentClass.h"
//#include <fstream>
//#include <math.h>

//const float s[5][5]={{1,0.5,0.5,0,0},{0.5,1,0.5,0,0},{0,0.5,1,0.5,0},{0.5,0,0.5,1,0},{0,0,0,0,0}};//打分矩阵
//const float s[5][5]={{1.0,-1,-2.0,-1,-2.0},{-1,1.0,-1,-2.0,-2.0},{-2.0,-1,1.0,-1,-2.0},{0,-2.0,0,1.0,-2.0},{-2.0,-2.0,-2.0,-2.0,0}};//打分矩阵
const float s[5][5]={{1,0,0,0,-2},{0,1,0,0,-2},{0,0,1,0,-2},{0,0,0,1,-2},{-2,-2,-2,-2,-2}};//打分矩阵
const float K = 1.23891;
const float L = 1.386294;

float tmpArray[PSS4LineLength][PSS4LineLength];
float ArrayMax(int a,int b,int& rangei,int& rangej);
float ShowMax(float a,float b,float c);
float S(char a,char b);

const float reducematrix=0;
const float gap=-1;


float ExtAlign(AlignmentResult &rar,int tb,int te,int pb,int pe)
{
	int i,j;
	float score=0;
	
	std::ofstream log("log1");
	log<<"Seed:"<<std::endl;
	for(i=rar.TBegin;i<rar.TEnd;i++)
		log<<rar.Text[i];
	log<<std::endl;
	for(i=rar.PBegin;i<rar.PEnd;i++)
		log<<rar.Pattern[i];
	log<<std::endl;
	log<<std::endl;
	for(i=tb;i<te;i++)
		log<<rar.Text[i];
	log<<std::endl;
	for(i=pb;i<pe;i++)
		log<<rar.Pattern[i];
	log<<std::endl;
	for(i=0;i<=(te-tb);i++)
	{
		for(j=0;j<=(te-tb);j++)
		{
			log<<tmpArray[i][j]<<"  ";
		}
		log<<std::endl;
	}
	

	for(i=0;i<=(te-tb);i++)
	{
		int j1=(i-(int)(reducematrix*(te-tb)))>0?(i-(int)(reducematrix*(te-tb))):0;
		int j2=(i+(int)(reducematrix*(te-tb))<(te-tb))?(i+(int)(reducematrix*(te-tb))):(te-tb);
		for(j=j1;j<=j2;j++)
		{
			if(i==0||j==0)
				tmpArray[i][j]=0;

			else
				tmpArray[i][j]=ShowMax(tmpArray[i-1][j-1]+S(rar.Text[tb+j-1],rar.Pattern[pb+i-1]),tmpArray[i][j-1]+gap,tmpArray[i-1][j]+gap);
		}
	}
	for(i=0;i<=(te-tb);i++)
	{
		for(j=0;j<=(te-tb);j++)
		{
			log<<tmpArray[i][j]<<"  ";
		}
		log<<std::endl;
	}
 	int rangei=pe-pb,rangej=te-tb;
	score=ArrayMax(pe-pb,te-tb,rangei,rangej);


	log.close();
	return score;
	
}

float AlignmentClass::FinalAlign(AlignmentResult &rar,int tb,int te,int pb,int pe)
{
	int i,j;
	float score=0;
	int tgapcount=0;
	int pgapcount=0;
	
	//int tmpArray[PSS4LineLength][PSS4LineLength];
	std::ofstream log("log1");
	
	for(i=0;i<=(pe-pb);i++)
	{
		//int j1=(i-(int)(reducematrix*(te-tb)))>0?(i-(int)(reducematrix*(te-tb))):0;
		//int j2=(i+(int)(reducematrix*(te-tb))<(te-tb))?(i+(int)(reducematrix*(te-tb))):(te-tb);
		for(j=0;j<=(te-tb);j++)
		{
			
			if(i==0||j==0)
				tmpArray[i][j]=0;
			else
				tmpArray[i][j]=ShowMax(tmpArray[i-1][j-1]+S(rar.Text[tb+j-1],rar.Pattern[pb+i-1]),tmpArray[i][j-1]+gap,tmpArray[i-1][j]+gap);
			  //tmpArray[i][j]=ShowMax(tmpArray[i-1][j-1]+S(rar.Text[tb+j-1],rar.Pattern[pb+i-1]),tmpArray[i][j-1]-gap,tmpArray[i-1][j]-gap);
		}
	}
	int rangei=pe-pb,rangej=te-tb;
	score=ArrayMax(pe-pb,te-tb,rangei,rangej);
    /*log<<"\t\t";
	for(j=0;j<(te-tb);j++)
		log<<rar.Text[j]<<"\t";
	log<<std::endl;*/
	/*for(i=0;i<=(pe-pb);i++)
	{
		if(i>0)
			log<<rar.Pattern[i-1]<<"\t";
		else
			log<<"\t";
		for(j=0;j<=(te-tb);j++)
		{
			log<<tmpArray[i][j]<<"\t";
		}
		log<<std::endl;
	}*/
 
	//for(i=(pe-pb),j=(te-tb);i>0&&j>0;)
	if(rangei<pe-pb)
	{
		//rar.Text=rar.gapinsert(rar.Text,te-tb,pe-pb-rangei,te-tb);
		for(int i=0;i<pe-pb-rangei;i++)
		{
			rar.TGapRecord.push_back(te-tb);
			tgapcount++;
		}
		//rar.Text[te-tb+pe-pb-rangei]='\0';
			
	}
	if(rangej<te-tb)
	{
		//rar.Pattern=rar.gapinsert(rar.Pattern,pe-pb,te-tb-rangej,pe-pb);
		for(int i=0;i<te-tb-rangej;i++)
		{
			rar.PGapRecord.push_back(pe-pb);
			pgapcount++;
		}

		//rar.Pattern[pe-pb+te-tb-rangej]='\0';
	}
	for(i=rangei,j=rangej;i>0&&j>0;)
	{
		float ttt=S(rar.Text[tb+j-1],rar.Pattern[pb+i-1]);
		if(tmpArray[i][j]==tmpArray[i-1][j-1]+ttt)
		{
			i--;
			j--;
		}
		else if(tmpArray[i][j]==tmpArray[i-1][j]+gap)
		{
			//rar.Text=rar.gapinsert(rar.Text,j-1,1,te-tb);
			rar.TGapRecord.push_back(j-1);
			tgapcount++;
			i--;
		}
		else if(tmpArray[i][j]==tmpArray[i][j-1]+gap)
		{
			//rar.Pattern=rar.gapinsert(rar.Pattern,i-1,1,pe-pb);			
			rar.PGapRecord.push_back(i-1);
			pgapcount++;
			j--;
			
		}
	}
	if(i>0)
	{
		//rar.Text=rar.gapinsert(rar.Text,0,i,te-tb);
		for(int tmpi=0;tmpi<i;tmpi++)
		{
			rar.TGapRecord.push_back(-1);
			tgapcount++;
		}

	}
	if(j>0)
	{
		//rar.Pattern=rar.gapinsert(rar.Pattern,0,j,pe-pb);
		for(int tmpj=0;tmpj<j;tmpj++)
		{
			rar.PGapRecord.push_back(-1);
			pgapcount++;
		}

	}
	//score=tmpArray[pe-pb][te-tb];
	
	rar.SimScore=score;
	rar.Similarity=rar.SimScore/float(te-tb+pgapcount+tgapcount);
	rar.TGapCount=tgapcount;
	rar.PGapCount=pgapcount;
	rar.gapinsert();
	return score;
}


void AlignmentClass::MainAlignment(AlignmentResult &rar)
{
	/*char *T=tar.Text;
	char *P=tar.Pattern;
	AlignmentResult rar;
	ResultAsignMent(&rar,&tar);*/
	bool RFlag=true,LFlag=true;
	int extcount=0;
	while(RFlag==true)//while(LFlag==TRUE||RFlag==TRUE)
	{
		if(LFlag==true)
		{
			//LeftExtend
			if(rar.IsPbegin()||rar.IsTbegin())
			{
				LFlag=false;
			}
			else
			{
				
				float Score=0;
				int Lext=rar.Lextern();
				
				Score=ExtAlign(rar,rar.TBegin-Lext,rar.TBegin,rar.PBegin-Lext,rar.PBegin);
			
				//int t1=(Score+rar.SimScore),t2=(Lext+(rar.PEnd-rar.PBegin)),t3=(rar.SimScore)+t1,t4=(rar.PEnd-rar.PBegin)+t2;
				
				if(Score/Lext<LocalSimDomain)
					LFlag=false;
				else
				{
					extcount++;
					//RFlag=true;
					rar.TBegin=rar.TBegin-Lext;
					rar.PBegin=rar.PBegin-Lext;
					//rar.SimScore=Score;
					//rar.Similarity=rar.SimScore/(rar.PEnd-rar.PBegin);
				}
			}
			if(RFlag==true)
			{
				
				//RExtend
				if(rar.IsPend()||rar.IsTend())
				{
					RFlag=false;
				}
				else
				{
					int Score=0;
					int Rext=rar.Rextern();

					Score=ExtAlign(rar,rar.TEnd,rar.TEnd+Rext,rar.PEnd,rar.PEnd+Rext);

					if(Score/Rext<LocalSimDomain)
						RFlag=false;
					else
					{
						extcount++;
						//LFlag=true;
						rar.TEnd=rar.TEnd+Rext;
						rar.PEnd=rar.PEnd+Rext;
						//rar.SimScore=Score;
						//rar.Similarity=rar.SimScore/(rar.PEnd-rar.PBegin);
					}					
				}
			}
		}
		
		else
		{
			if(RFlag==true)
			{
				
			//	RExtend();
				if(rar.IsPend()||rar.IsTend())
				{
					RFlag=false;
				}
				else
				{
					int Score=0;
					int Rext=rar.Rextern();

					Score=ExtAlign(rar,rar.TEnd,rar.TEnd+Rext,rar.PEnd,rar.PEnd+Rext);

					if(Score/Rext<LocalSimDomain)
						RFlag=false;
					else
					{
						extcount++;
						//LFlag=true;
						rar.TEnd=rar.TEnd+Rext;
						rar.PEnd=rar.PEnd+Rext;
						//rar.SimScore=Score;
						//rar.Similarity=rar.SimScore/(rar.PEnd-rar.PBegin);
					}
				}				
			}
		}
	}
	if(extcount>=3)
	{
		rar.SimScore=FinalAlign(rar,rar.TBegin,rar.TEnd,rar.PBegin,rar.PEnd);
		rar.Similarity=rar.SimScore/(rar.PEnd-rar.PBegin);
		float m=rar.PatternCount;
		float n=rar.TextCount;
		float x=rar.SimScore;
		rar.p_value=1-exp(-K*m*n*exp(-L*x));
	}
	else
	{
		rar.SimScore=-1;
		rar.Similarity=-1;
		rar.p_value=-1;
	}
	//return rar;
}

int AlignmentResult::Lextern()
{
	int i=(TBegin-ExtendStep)>0?ExtendStep:TBegin;
	int j=(PBegin-ExtendStep)>0?ExtendStep:PBegin;
	i=(i<=j?i:j);
	return i;
}
int AlignmentResult::Rextern()
{
	int i=(TEnd+ExtendStep)<TextCount?ExtendStep:TextCount-TEnd;
	int j=(PEnd+ExtendStep)<PatternCount?ExtendStep:PatternCount-PEnd;
	i=(i<=j?i:j);
	return i;
}

bool AlignmentResult::IsTbegin()
{
	if(TBegin==0)
		return true;
	else
		return	false;

}

bool AlignmentResult::IsTend()
{
	if(TEnd==TextCount)
		return true;
	else
		return	false;

}
bool AlignmentResult::IsPend()
{
	if(PEnd==PatternCount)
		return true;
	else
		return	false;

}
bool AlignmentResult::IsPbegin()
{
	if(PBegin==0)
		return true;
	else
		return	false;

}

void AlignmentClass::ResultAsignMent(AlignmentResult& rar,AlignmentResult& tar)
{
	rar.PBegin=tar.PBegin;
	rar.PEnd=tar.PEnd;
	rar.TBegin=tar.TBegin;
	rar.TEnd=tar.TEnd;

	rar.Similarity=tar.Similarity;
	rar.SimScore=tar.SimScore;

	rar.PatternCount=tar.PatternCount;
	rar.TextCount=tar.TextCount;
	rar.Text=tar.Text;
	rar.Pattern=tar.Pattern;
}
float S(char a,char b)//打分矩阵
{
	int a1,b1;
	switch(a)
	{
		case '0':
			 a1=0;break;
		case '1':
			 a1=1;break;
		case '2':
			 a1=2;break;
		case '3':
			 a1=3;break;
		case ' ':
			 a1=4;break;
	}
	switch(b)
	{
		case '0':
			 b1=0;break;
		case '1':
			 b1=1;break;
		case '2':
			 b1=2;break;
		case '3':
			 b1=3;break;
		case ' ':
			b1=4;break;
	}
	if(a1==-1||b1==-1)
		return 0;
	float t=s[a1][b1];
	return t;
	
}
float ShowMax(float a,float b,float c)
{
	float tmp=a>b?a:b;
	tmp=tmp>c?tmp:c;
	return tmp;
}

void AlignmentClass::test(AlignmentResult &tmpAR)
{
	std::ofstream OutFile;
	OutFile.open("odata\\test.out",std::ios::out|std::ios::app);
	//int tmpExt=Result.PEnd-Result.PBegin;
	int i,j;
	for(i=tmpAR.PBegin;i<tmpAR.PEnd;i++)
	{
		for(j=tmpAR.TBegin;j<tmpAR.TEnd;j++)
		{
			if(tmpArray[i][j]>=0)
				OutFile<<tmpArray[i][j]<<" ";
			else
				OutFile<<"*"<<" ";
		}
		OutFile<<std::endl;
	}
	

}
float ArrayMax(int a,int b,int &rangei,int &rangej)
{
	float max;
	if(a<=0||b<=0)
	{
		rangei=0;
		rangej=0;
		return max=tmpArray[0][0];
	}
	else
	{
		max=tmpArray[a][a];
		rangei=a;
		rangej=b;
		for(int i=1;i<a;i++)
		{
			if(max<tmpArray[a-i][b])
			{
				rangei=a-i;
				max=tmpArray[a-i][b];
			}
		}
		for(int j=1;j<b;j++)
		{
			if(max<tmpArray[a][b-j])
			{
				rangej=b-j;
				max=tmpArray[a][b-j];
			}
		}


		return max;
	}
	
}

void AlignmentResult::gapinsert()
{


	int PosMove=TGapRecord.size();
	if(!TGapRecord.empty())
	{
		int tmp=TGapRecord[PosMove-1];
		while(TGapRecord[PosMove-1]==-1)
		{
			AlignedT.push_back('-');
			if(PosMove>0)
				PosMove--;
			else
				break;

		}
	}
	for(int i=0;i<TextCount;i++)
	{
		AlignedT.push_back(Text[i]);
		int tmp=TGapRecord[PosMove-1];
		if(!TGapRecord.empty())
		{
			if(PosMove>0)
			{

				while(TGapRecord.at(PosMove-1)==i)
				{
					AlignedT.push_back('-');
					if(PosMove>0)
						PosMove--;
					else
						break;
				}
			}
		}
	}
	for(i=0;i<PosMove;i++)
		AlignedT.push_back('-');

	PosMove=PGapRecord.size();
	if(!PGapRecord.empty())
	{
		while(PGapRecord[PosMove-1]==-1)
		{
			AlignedP.push_back('-');
			if(PosMove>0)
				PosMove--;
			else
				break;

		}
	}
	for(i=0;i<PatternCount;i++)
	{
		AlignedP.push_back(Pattern[i]);
		if(!PGapRecord.empty())
		{

			if(PosMove>0)
			{
				while(PGapRecord.at(PosMove-1)==i)
				{
					AlignedP.push_back('-');
					if(PosMove>0)
						PosMove--;
					else
						break;
				}
			}
		}
	}
	for(i=0;i<PosMove;i++)
		AlignedP.push_back('-');


}

⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -