main.cpp

Pairwise Sequence alignment Algorithm Simulation. Read 2 input strings and then calculate the Pairw

#include<iostream>
#include<cstdio>
#include<iomanip>
#include <conio.h>



#include <io.h>
#include <fcntl.h>
#include <stdlib.h>


using namespace std;




	void traceback(int i, int j, int **b,char x[], char y[])
	{
		if(i==0 && j==0)			//Finish trace back
			return;

	
		if(b[i][j]==1)				//The value is from diagonal
		{
			traceback(i-1, j-1, b, x, y);
			cout<<" "<< x[i-1]<<"       "<<y[j-1]<<endl;

			FILE *fp1,*fp2;
			fopen_s(&fp1,"xprime.txt","a");	//Open x'
			fopen_s(&fp2,"yprime.txt","a");	//Open y'

			fprintf(fp1,"%c",x[i-1]);		//Save result to x'
			fprintf(fp2,"%c",y[j-1]);		//Save result to y'

			fclose(fp1); 
			fclose(fp2); 
		}
		else
			if(b[i][j]==2)			//The value is from top
			{
				traceback(i-1, j, b,x, y);
				cout<<" "<< x[i-1]<<"       "<<"-"<<endl;

				FILE *fp1,*fp2;
				fopen_s(&fp1,"xprime.txt","a");
				fopen_s(&fp2,"yprime.txt","a");

				fprintf(fp1,"%c",x[i-1]);
				fprintf(fp2,"%c",'-');

				fclose(fp1); 
				fclose(fp2); 
			}
			else					//The value is from left
			{
				traceback( i, j-1,b, x, y);
				cout<<" "<< "-"<<"       "<<y[j-1]<<endl;

				FILE *fp1,*fp2;
				fopen_s(&fp1,"xprime.txt","a");
				fopen_s(&fp2,"yprime.txt","a");

				fprintf(fp1,"%c",'-');
				fprintf(fp2,"%c",y[j-1]);

				fclose(fp1); 
				fclose(fp2); 
			}
	
	}


int  main()
{

	int lenX;							// Length of X
	int lenY;							// Length of Y
	int lenResult=0;					// Length of X' and Y'




	cout<<"Welcome to use this program!"<<endl<<endl;
	cout<<"This program will read two input sequences from files, calculate the C table and B table, save the results to files and then display the results X', Y'."<<endl<<endl;

	cout<<"Press any key to continue..."<<endl<<endl;
	_getch();
	


/////////////////////////////////////////////
//
// Read and store the file x.txt and y.txt
//
/////////////////////////////////////////////


	// Get the length of X
	int fh;
	_sopen_s( &fh, "x.txt", _O_RDONLY, _SH_DENYNO, 0 );

	lenX = _lseek( fh, 0L, SEEK_END );
	if( lenX == -1L )
		perror( "_lseek to end failed" );
	else
		cout<<"The length of X is "<<lenX<<endl;

	_close(fh);	


/////////////////////////////////////////////////////////////////////////

	// Open x.txt and store the input to x[]
	FILE *fp;
	fopen_s(&fp,"x.txt","r");

	char *x;
	x = new char[lenX+1];		// Allocate space for array X

 	fscanf_s(fp,"%s",x);		// Store the input file to x[]
	fclose(fp);					//close file


/////////////////////////////////////////////////////////////////////////


	// Get the length of Y

	_sopen_s( &fh, "y.txt", _O_RDONLY, _SH_DENYNO, 0 );

	lenY = _lseek( fh, 0L, SEEK_END );
	if( lenY == -1L )
		perror( "_lseek to end failed" );
	else
		cout<<"The length of Y is "<<lenY<<endl<<endl;

	_close(fh);	

/////////////////////////////////////////////////////////////////////////

	// Open y.txt and store the input to y[]

	fopen_s(&fp,"y.txt","r");

	char *y;
	y = new char[lenY+1];		// Allocate space for array Y

 	fscanf_s(fp,"%s",y);		// Store the input file to y[]
	fclose(fp); //close file




/////////////////////////////////////////////
// Print out the input X and Y
/////////////////////////////////////////////


	cout<<"X = ";
	for(int i=0; i < lenX; i++) // Print X
	{
		cout <<x[i]<<" ";
	}
	cout<<endl;
	
	cout<<"Y = ";
	for(int j=0; j < lenY; j++) // Print Y
	{
		cout <<y[j]<<" ";
	}
	cout<<endl;



/////////////////////////////////////////////
// Create c table and b table
/////////////////////////////////////////////


	int **c = new int*[lenX+1];        // Allocate space for array c[][]
	c[0] = new int[(lenX+1)*(lenY+1)];
	for(int i=0;i<(lenX+1);i++)
	{
		c[i]=c[0]+i*(lenY+1);
	}


	int **b = new int*[lenX+1];        // Allocate space for array b[][]
	b[0] = new int[(lenX+1)*(lenY+1)];
	for(int i=0;i<(lenX+1);i++)
	{
		b[i]=b[0]+i*(lenY+1);
	}


	for(int i=0; i <= lenX; i++) 		// Initialize c table. first colum = -2
	{
		c[i][0] = -2;
	}


	for(int j=0; j <= lenY; j++)		// Initialize c table. first row = -2
	{
		c[0][j] = -2;
	}


	for(int i=0; i <= lenX; i++)		 // Initialize b table. first colum = 2
	{
		b[i][0] = 2;
	}

	for(int j=0; j <= lenY; j++)		 // Initialize b table. first row = 3
	{
		b[0][j] = 3;
	}

	
	// Calculate c table and b table

	for(int i=1; i <= lenX; i++)
		for(int j=1; j <= lenY; j++)
		{
			int dia,top,lef;

			//Calculate the value of diagonal, top and left

			if(x[i-1]==y[j-1])		//Calculate the value of diagonal
			{	
				dia=c[i-1][j-1]+1;
			}
			else
				dia=c[i-1][j-1]-1;

			top=c[i-1][j]-2;		//Calculate the value of top
			lef=c[i][j-1]-2;		//Calculate the value of left

			// Find the max{ dia, top, lef }
			if((dia >= top) && (dia >=lef))
			{
				c[i][j]=dia;		//If the value of c[i][j] is from diagonal
				b[i][j]=1;			//Set b[i][j] = 1
			}
			else
				if(top >= lef)
				{
					c[i][j]=top;	//If the value of c[i][j] is from top
					b[i][j]=2;		//Set b[i][j] = 2
				}
				else
				{
					c[i][j]=lef;	//If the value of c[i][j] is from left
					b[i][j]=3;		//Set b[i][j] = 3
				}
		}

/*
	//Print out C table
	cout <<endl<<"C Table:"<<endl;
	for(int i=0; i <= lenX; i++)
	{
		for(int j=0; j <= lenY; j++)
		{
			cout << setw(3) << c[i][j];
		}
		cout<<endl;
	}



	//Print B table
	cout <<endl<<"B Table:"<<endl;
	for(int i=1; i <= lenX; i++)
	{
		for(int j=1; j <= lenY; j++)
		{
			if(b[i][j] ==1)			//Diagonal
				cout << "\\  ";		
			else
				if(b[i][j] ==2)		//Top
					cout<<"|  ";
				else
					cout<<"-  ";	//Left
		}
		cout<<endl;
	}

*/

/////////////////////////////////////////////
//
// Trace back
//
/////////////////////////////////////////////



	cout<<endl<<endl<<" X' <->  Y'"<<endl;
	cout<<"          "<<endl;

	traceback(lenX,lenY,b,x,y);





	cout<<"\n\n\nPress any key to exit this program.\n";

	_getch();

	delete x;
	delete y;

	return 0;

}