/*

This is a program to create genome Sequences and align genome subsequences using Longest Common Subsequence (LCS).
The program runs on a Unix platform and can be run by entering the size of the sequence and the number of 
subsequences form the command prompt.

@Author : Sara Nasser

@Date Created: January 26, 2006

@Date Modified: February 28, 2006

@Version :  3.2 

@Modified by: Sara Nasser

(c) This code is created by the bioinformatics research group in the Department of 
Computer Science at the University of Nevada Reno.
 
*/

//System includes
#include<fstream.h>
#include <iostream.h>
#include <iostream.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>

//user directives
#define display 0	//Display turned on when set to 1 and off when set to 0
//#define NSQ  100  // no of sub sequences


/* Function prototypes */

void Subsequence(int, int);
void splicer();
void LCS(char *, char*, char *, int, int&); 
void substring(char* ,char* , int , int );
void mergeSub(char* ,char* ,char* , char * , int);
void align(char**, int, char*, ofstream&);
void display_html_original(ofstream& , char **, int);
void display_html_1(ofstream&, char** );
void display_html_2(ofstream&, char* ,char*, float);
void display_html_merged(ofstream&, char *);
void display_error(ofstream& , int , int);


int SIZE=300; // the size of the genome sequence

/*****************************************************************

The main functions reads two variables from the command prompt, 

  @return: int, 0 if the programs runs successfully, 1 if it does not

  @param :argc the number of arguments, argv the arguments as strings

******************************************************************/
int main(int argc, char* argv[])
{	
	
	char newSeq[SIZE+1], nextSeq[SIZE+1], overlap[SIZE+1], mergedSeq[SIZE+1];
	char mainSeq[SIZE+1];

	int noSeq;
	int overlapSize,oldOverSize=0;
	int type;

	//int NSQ=100;
	
	SIZE=atoi(argv[1]);  // get the size of the main sequence from the user

	int NSQ=atoi(argv[2]); //get the number of subsequences desired

	//if sequence size is small and number of subsequences are less than two display an error message
	if(SIZE < 10 || NSQ <2 )
	{
	if(display ==1)
	{
 //If user chooses display, a html file is created that contains the results.
        ofstream fout2;
	 fout2.open("genome.html");
	 display_error(fout2, SIZE, NSQ);	
	}		
	return 0;
	}

	cout<<SIZE<<"  "<<NSQ<<" "; 
	
	time_t seconds1, seconds2;
	float totaltime;

	int i,j;
	ifstream ifp;
	ifp.open("subseq.dat"); // a file containing the sequences and subsequences which will be used for alignment, this feature will be used later!

	//If user chooses display, a html file is created that contains the results.
	ofstream fout2;
	if(display==1)
	  {
	 //If user chooses display, a html file is created that contains the results.
        //ofstream fout2;


	  // getenv("HOME");
	  //fout2.open("/public_html/code/genome.html");
	    fout2.open("genome.html");	  
	}
	//create data structure to store the sub sequences
	char **data;
	data= new char*[NSQ]; //allocate memory to data NSQxSIZE  
	for(i=0; i<NSQ; i++)
	{
		data[i]= new char[SIZE+1];
	}


	// generate test data
	Subsequence(SIZE,NSQ);

	//start timer
	seconds1 = time (NULL);

	//read the main sequence from the file
	ifp.getline(mainSeq, SIZE+1, '\n');

	//read all thesubsequences into the variable data
	for( i=0; i<NSQ; i++)
	{
		ifp.getline(data[i], SIZE+1, '\n');
	}

	//If display is turned on, the display file is initialized
	if(display==1)
	   display_html_1(fout2, data);	
	

	//aling sequences
	align(data, NSQ, mergedSeq, fout2);

	//stop timer
	seconds2 = time (NULL);

//	cout<<"\n\n ------------Main Sequence---------------\n "<<mainSeq<<endl;
//	cout<<"\n\n ------------Final Sequence--------------------\n "<<mergedSeq;

//	cout<<"\n\n ------------Length of final sequence: "<<strlen(mergedSeq);

	cout<<strlen(mergedSeq);

	//check whether the aligned sequences matches the main sequence and give message to user
	if(!strcmp(mainSeq,mergedSeq))
	{
		//cout<<"\n\nThey match!!\n";
		cout<<" ";
	}
	else
	{
		//cout<<"\n\nCould not get the entire sequence but the largest!!\n";
		cout<<" ";
	}

	totaltime=(seconds2-seconds1); //   /(60.0);
	//cout<<"\nTotal time for assembly is "<<(seconds2-seconds1)<<" seconds "<<endl;
	cout<<totaltime<<endl;

	//display the final results
	if(display==1)
	   display_html_2(fout2, mainSeq, mergedSeq, totaltime);
	
	//cout<<"1="<<NSQ;

		
	
//free space
/*
for(i=0; i<=NSQ; i++)
{  
	delete[] data[i];
//	delete[] mergeRes[i];
}
delete[] data;
//delete[] mergeRes;

*/
return 0;

}

/**********************************************

Function to chop the sequence into subsequences

   @return: void

   @param : seqSize is the length of the sequence
			noSubSeq the number of subsquences
*************************************************/

void Subsequence(int seqSize, int noSubSeq)
{
	
	srand(time(NULL));  //set the random seed

	char nucleo[4]={'A', 'T', 'C', 'G'};	//array for the 4 nucleotides
	char seq[SIZE+1];                       //array that has the actual sequences 
	int start, endse, temp;
	
	//output data to a file
	ofstream fout;
	fout.open("subseq.dat");				//file to store the subsequences

//	cout<<"Sequence .........................\n";

	int i;
	for(i=0; i<seqSize; i++)				//generate the main sequence
	{
		seq[i]=nucleo[rand()%4];	
		//cout<<seq[i];
		fout<<seq[i];
	}
	
	fout<<endl;

	//fout<<noSubSeq<<endl;

	for(i=0; i<noSubSeq; i++)		//generate the subsequences
	{
//		cout<<"\nSub Sequence "<<i<<" .........................\n";
		start=rand()%seqSize; 
		endse=rand()%seqSize;
		if(start > endse)
		{
			temp=start;
			start=endse;
			endse=temp;
		}
		while(start==endse)// regenerate sub sequence if empty
		{
			start=rand()%seqSize; 
			endse=rand()%seqSize;
			
			if(start > endse)
			{
				temp=start;
				start=endse;
				endse=temp;
			}

		}
//		cout<<"Start "<<start<<" End "<<endse<<endl;

		int j;
		for( j=start; j<=endse; j++)	//store the subsequence in the file
		{ 
	//		cout<<seq[j];
			fout<<seq[j];
		}
		fout<<endl;
	}
	
	fout.close();

}

/**********************************************

Function to aling sequences, exact alignment 
 
   @return: void

   @param :data character array with the subsequences
			nsq the number of subsquences
			mRes  character array for final merged sequence			
			fout2 the output file stream for display

**********************************************/	
void align(char **data, int nsq, char *mRes, ofstream& fout2)
{
		
	char newSeq[SIZE+1], nextSeq[SIZE+1], overlap[SIZE+1], mergedSeq[SIZE+1],maxSeq[SIZE+1];
	int noSeq;
	int overlapSize,oldOverSize=0;
	int type;
	int i,j, maxlen=0;


	
	mRes[0]='\0';
	i=0;
/*		
	while(i<nsq)
	{
	  maxSeq[0]='\0'; 
		
		//the original sub sequence
	  //	memcpy(maxSeq,data[i],strlen(data[i])*sizeof(char)+1);
		
		
		//next one onwards
		for( j=i+1; j<nsq; j++)
		{
			//find LCS
			LCS(data[i],data[j],overlap,2, type);
			
			overlapSize=strlen(overlap);
			
			if(overlapSize>0 && type>0)	//has a match, creates an empty mergedSubst if not checked, to be used later to define size of overlap  
			{
				mergeSub(data[i], data[j],overlap, mergedSeq, type);
			
				//show only new sequences
				if((strlen(mergedSeq) > strlen(data[i])) && (strlen(mergedSeq) > strlen(data[j])))
				{
					//pick up the longest merged subsequence   
					if(strlen(mergedSeq) >= strlen(maxSeq))
					{
						memcpy(maxSeq,mergedSeq,strlen(mergedSeq)*sizeof(char)+1);
					}
				}
	
			}
		
	
		
		}//for

		//copy the longest sub sequence
		memcpy(data[i], maxSeq, strlen(maxSeq)*sizeof(char)+1);

		i++;
	}//while
*/
	
	////Combine the results into one final sequence
	j=0;
	int oldmax;
	maxlen=0;
	while((j < nsq) && (maxlen<SIZE))
	{
	//  cout<<"in while....";
	        oldmax=maxlen;
		maxSeq[0]='\0';
			
		 //memcpy(maxSeq,data[i],strlen(data[i])*sizeof(char)+1);

		if(strlen(data[j])!=0)
		{
		  for(i=0; i<nsq; i++)
		    {
		      if(strlen(data[i])!=0)
			{
				//find LCS
			  LCS(data[j],data[i],overlap,2, type);
			  
			  overlapSize=strlen(overlap);
			  
			  if(overlapSize>0 && type>0)	//has a match, creates an empty mergedSubst if not checked, to be used later to define size of overlap  
			    {
			      mergeSub(data[j], data[i],overlap, mergedSeq, type);
			      
			      //show only new sequences
			      if((strlen(mergedSeq) > strlen(data[i])) && (strlen(mergedSeq) > strlen(data[j])))
				{
				  
				  //pick up the longest merged subsequence   
				  if(strlen(mergedSeq) >= strlen(maxSeq))
				    {
				      //if(display)
				      //display_html_original(fout2,data, nsq);
				      
				      memcpy(maxSeq,mergedSeq,strlen(mergedSeq)*sizeof(char)+1);
				      
				      if(strlen(mergedSeq) >= maxlen)
					maxlen=strlen(maxSeq);
				    }
				}//if
			      
			    }//if
			}//if	
		      
		    }//for
		  
		  //copy the longest sub sequence

		  //If subsequence does not match any of the other merged/unmerged
		  //subsequences we make it null, so that we don't use it for
		  //any further alignment.
		//  if(strlen(maxSeq)==0)
		 //   data[j]=NULL;
		 // else
		
		    data[j]=maxSeq;

		  //cout<<"new res="<<mergeRes[j];
		  
		  
		  
		} //if
		
		
		if(strlen(data[j]) >= strlen(mRes))
		  {
		    if(display==1)
		      display_html_original(fout2, data, nsq);
		    
		    memcpy(mRes,data[j],strlen(data[j])*sizeof(char)+1);
		    //cout<<"Max so far"<<mRes<<endl;
		    maxlen=strlen(mRes);
		    //cout<<"Max len"<<maxlen;
		    
		    if(display==1)
		      display_html_merged(fout2, mRes);
		  }
		j++;
		
		//if the string is changing we continue
		if((j==nsq) && (maxlen > oldmax))
		  {
		    j=0;
		  }
		
	}//while
	
}

/**********************************************************************************

Function to put back the subsequences, using LCS -Dynamic programming
Conventions used -1 for up, -2 for left and -3 for up and left in the table

  @return: void

  @param :x and y are the sub sequences to be compared,
		  z is the overlap sub sequences, 
		  style is the type of overlap 2 :contiguous and 1:non-contiguous.
		  type is a return variable 1 if z is a suffix of x prefix of y
									2 if z is a suffix of y prefix of x
									0 if z is not a suffix or prefix of x,y
***********************************************************************************/
void LCS(char *x, char *y, char *z, int style, int &type)
{

//variables	
int i,j, xstart, ystart, xend, yend;

int x_l, y_l;


x_l=strlen(x);   //lenght of string x and y
y_l=strlen(y);

//if strings are not empty
if(x_l!=0 && y_l!=0)
{

int **c, **b;		//allocate space
c= new int* [x_l+1];
b= new int* [x_l+1];

for(i=0; i<=x_l; i++)
{
	b[i]=new int[y_l+1];
	c[i]=new int[y_l+1];

	for(j=0; j<=y_l; j++)
		b[i][j]=0;
}


for(i=0; i<=x_l; i++)   //initialize the tables b and c
{
		c[i][0]=0;
		b[i][0]=-1;   //-1 indicates up
}


for(i=0; i<=y_l; i++)   //initialize the tables b and c
{
		c[0][i]=0;
		b[0][i]=-2;  //-2 indicates left
}


//dynamic programming, constructs the table that can be used to reconstruct the LCS

for(i=1; i<=x_l; i++)
	for(j=1; j<=y_l; j++)
	{
		if(x[i-1]==y[j-1])  //two characters in the strings match
		{
			c[i][j]=c[i-1][j-1]+1;
		
				b[i][j]=-3;           //indicates up and left
		
		}
		else if(c[i-1][j]>=c[i][j-1])
		{
			c[i][j]=c[i-1][j];
			b[i][j]=-1;			//up
		}
		else
		{
			c[i][j]=c[i][j-1];
			b[i][j]=-2;			//left
		}
	
	}

	//the largest no is the length of the common string
	int len=c[x_l][y_l];
	int start, stop=0, max, oldstart, oldstop, ctr, len2=0;

	int jcopy;

	
	z[len]='\0';

	i=x_l; j=y_l;

//back track
	oldstop=0; stop=0; ctr=0;
	oldstart=0; yend=0;
	start=i; jcopy=j;



while(i>0 || j>0)
{
	if(b[i][j]==-3)
	{
		i--; j--;
		len--;
		z[len]=x[i]; //copy the character

		
		if(ctr==0)		//save the starting position
		{
			start=i;
			jcopy=j;
			stop=0;
			ctr=1;
		}
		else
		{

			stop++;

			if(i==0)
			{
				if(stop>=oldstop)
				{
					oldstop=stop;
					oldstart=start;
					yend=jcopy;
				}

			}
		}
	}//if -3
	else if(b[i][j]==-1)
	{
		i--;
	
		if(stop>=oldstop)
		{
			oldstop=stop;
			oldstart=start;
			yend=jcopy;

		} 
		ctr=0;
		
	}//else if -1
	else if(b[i][j]==-2)
	{
		j--;

		if(stop>=oldstop)
		{
			oldstop=stop;
			oldstart=start;
			yend=jcopy;
		}	
		ctr=0;
		
	}//else if -2

}//while



int i1,j1, count=0, oldcount=0;

i=oldstart-oldstop+1;
j=yend-oldstop+1;

oldstart=xstart=i;
ystart=j;


//to find the LCS from the starting position
if(style==2)
{
	while( i<=x_l && j<=y_l &&(b[i][j]==-3))
	{

		i++;
		j++;
	count++;
	
	}
	xend=oldstart=i-1;
	yend=j-1;
	oldstop=count;
}



if(oldstop ==0)
{
	z[0]='\0';
	type=0;   // no match
}

else if(oldstop>0 && style==2) //if there is a match
{	
	xstart=oldstart-oldstop;
	xend=oldstart;

	ystart=yend-oldstop;
	
	oldstart--;
	substring(z, x, xstart, oldstart); //we only want it if its more than 1 char
	
	//Determine if the common substring is at the ends
	if((xend ==x_l) &&  (ystart==0)) //suffix of x prefix of y
	{

		type = 1;
	}
	else if((xstart==0) && (yend==y_l))  //suffix of y prefix of x
	{

		type = 2;
	}
	else
		type = 0; // not a suffix or prefix
}




//free space

for(i=0; i<=x_l; i++)
{  
	delete[] b[i];
	delete[] c[i];
}
delete[] b;
delete[] c;

}// if x_l !=0 ....

  }

/************************************************************************

Function to merge two subsequences based on 
an exact match at the edges.

   @return: void

  @param :	x,y strings to be merged
			overlap string LCS	
			z result merged string
			type 1 if overlap is a suffix of x prefix of y
				 2 if overlap is a suffix of y prefix of x
				 0 if overlap is not a suffix or prefix of x,y

**************************************************************************/
void mergeSub(char* x,char* y,char* overlap, char * z, int type)
{
	int i,j;
	int len_m=strlen(overlap);
	int len_x=strlen(x);
	int len_y=strlen(y);

	if(type==1) // suffix of x prefix of y
	{
		strcpy(z,x); //copy x to new string
		for(i=len_x, j=len_m; j<len_y; i++,j++)
		{
			z[i]=y[j];
		}
		z[i]='\0';

	}
	else if (type ==2)  //suffix of y prefix of x
	{
		strcpy(z,y); //copy x to new string
		for(i=len_y, j=len_m; j<len_x; i++,j++)
		{
			z[i]=x[j];
		}
		z[i]='\0';
	}
	else if (type ==0)  //empty string
	{
		z[0]='\0';
	}

}

/**********************************************

Function to put back the subsequences-Not in Use

**********************************************/
void splicer(int type)
{


}
/***********************************************************************
Function to copy a string at the end of another

 @return: void

  @param :x,y strings ot be merged
		  start position in y
		  len the characters to be copied 	
**********************************************************************/
void substring(char* y,char* x, int start, int len)
{
  int i,j;
	for( i=start, j=0; i<=len; i++,j++)
	{
		y[j]=x[i]; 
	}
	y[j]='\0';
}

/*******************************************************

  This function creates a html display, This is one of the  4
  functions for display. It initializes the html file

	@return: void

	@param : fout2 is the ofstream instance which points to the html file,
			 data array with subsequences	
********************************************************/


void display_html_1(ofstream& fout2, char **data)
{

	int i;
	fout2<<"<html><head></head><body>";


	//display_html_original(fout2,data);


}
/*******************************************************

	This function creates a html display, This is one of 
	functions for display. It finalizes the html file

	@return: void

	@param :fout2 is the ofstream instance which points to the html file,
			orig the original main sequence
			merge the final merged sequence
			time in seconds, total time for alignment

********************************************************/
void display_html_2(ofstream& fout2, char *orig, char *merge, float time)
{	
	//also display time, length etc....
	fout2<<"<hr><br><h3>";

	fout2<<"<br>Original Genome Sequence:"<<"<u>"<<orig<<"</u>";
        fout2<<"<br>Length of Sequence: <u>"<<strlen(orig)<<"</u>";

	fout2<<"<br>Final Genome Sequence&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;:"<<"<u>"<<merge<<"</u>";
	fout2<<"<br>Length of Sequence: <u>"<<strlen(merge)<<"</u>";
	fout2<<"<br>Time:"<<"<u>"<<time<<"</u> seconds";
	fout2<<"</body></html>";
}

/*******************************************************
*
	This function creates a html display, This is one of the  5
	functions for display. It writes the sequences to the html file

	@return: void

	@param :
			fout2 is the ofstream instance which points to the html file,
			data array with subsequences
			NSQ the number of subsquences
********************************************************/

void display_html_original(ofstream& fout2, char **data, int NSQ)
{
	int i,j;
	int len;
//	char seq[DATA+1];

	fout2<<"<hr><h6>";
	fout2<<"<table><tr>";

	for(i=0; i<NSQ; i++)
	{
		fout2<<"<td>";
		if(strlen(data[i])<16 && strlen(data[i])> 0)
		{
			
			if(i%2)
			{
				fout2<<"<table border=2 bordercolor=\"black\"><tr><td>";
				fout2<<data[i]<<"</td></tr></table>";
			}
			else
			{
				fout2<<"<table border=2 bordercolor=\"red\"><tr><td>";
				fout2<<data[i]<<"</td></tr></table>";
			}
		}
		else if(strlen(data[i])>=16)
		{
			len=strlen(data[i]);
			if(i%2)
			{
				
				fout2<<"<table border=2 bordercolor=\"black\"><tr><td>";
				fout2<<data[i][0]<<data[i][1]<<data[i][2]<<data[i][3];
				fout2<<data[i][4]<<data[i][5]<<data[i][6]<<data[i][7];
				fout2<<data[i][8]<<"..."<<data[i][len-8]<<data[i][len-7]<<data[i][len-6];
				fout2<<data[i][len-4]<<data[i][len-3]<<data[i][len-2]<<data[i][len-1]<<"</td></tr></table>";
		//		fout2<<"<br><br>-------"<<data[i]<<"-----";
			}	
			else
			{
				fout2<<"<table border=2 bordercolor=\"red\"><tr><td>";
				fout2<<data[i][0]<<data[i][1]<<data[i][2]<<data[i][3];
				fout2<<data[i][4]<<data[i][5]<<data[i][6]<<data[i][7];
				fout2<<data[i][8]<<"..."<<data[i][len-8]<<data[i][len-7]<<data[i][len-6];
				fout2<<data[i][len-4]<<data[i][len-3]<<data[i][len-2]<<data[i][len-1]<<"</td></tr></table>";
				//fout2<<data[i]<<"</td></tr></table>";
			
			}	
		}
		fout2<<"</td>";
		//	fout2<<"&nbsp\;&nbsp\;&nbsp\;";
	}
	fout2<<"</tr></table>";
}

/*****************************************************************************

	This function creates a html display, This is one of the  5
	functions for display. It writes the final string to the html file

	@return: void

	@param :fout2 is the ofstream instance which points to the html file
			merge the final merged sequence
	
******************************************************************************/
void display_html_merged(ofstream& fout2, char *merged)
{
	
	int len;
	len=strlen(merged);

	fout2<<"<hr><h4>";
	fout2<<"<br>";
	fout2<<"<table><tr><td color=blue>";
	fout2<<"The longest Subsequence found so far...<br>"<<merged<<"<br>"<<"Lenght ="<<len<<"<br>";

}

/************************************************************************
*   This function displays an error in the browser if the size is small or
*   the number of subsequences areless
* 
    @return: void

	@param :fout2 is the ofstream instance which points to the html file,
			SIZE is the length of the sequence
			NSQ the number of subsquences
*
********************************************************/

void display_error(ofstream& fout2, int SIZE, int NSQ)
{
 //also display time, length etc....
        fout2<<"<html><head></head><body>";

	fout2<<"<hr><br><h3>";

        fout2<<"<br>The Program did not run"<<"<u>"<<"The reason is below"<<"</u>";
        fout2<<"<br>Length of Original Sequence: <u>"<<SIZE<<"</u>";

        fout2<<"<br>Number of SubSequences: <u>"<<NSQ<<"</u>";
        fout2<<" Length of Original Sequence should be greater than 10 and number of subsequences should be more than 2";
        fout2<<"</body></html>";

}