/code_comm/checkpcity.cpp

// checkpcity.cpp : Defines the entry point for the console application.

//



//#include "stdafx.h"

#include<stdio.h>

#include<fstream>

#include<iostream.h>

#include<math.h>

#include<string.h>

#include<stdlib.h>

//#include<comdef.h>

#include<ctype.h>



//#include "WindowAnalysis.h"

//#include "PowerSpectrum.h"

#include<vector>

#include<complex>

using namespace std;

class PowerSpectrum;







//Included the Powerspectrum module







class PowerSpectrum  

{

public:

	PowerSpectrum();

	virtual ~PowerSpectrum();

	char * strCRegion;

	double *s;

	double *x;

	//temp check

	double *x1;

	double *xp;

    double *xp1;

	double *xr;

	double *xr1;



	int  doPowerSpectrumAnalysis(); 

	void	setSequence(const char * seq);

	

};







PowerSpectrum::PowerSpectrum()

{

	s=NULL;

	x=NULL;

	strCRegion=NULL;



}



PowerSpectrum::~PowerSpectrum()

{

	delete []s;

	delete []x;

    delete []strCRegion;



}



void PowerSpectrum::setSequence(const char * seq)

{



delete []strCRegion;



strCRegion=new char[ strlen(seq)+1 ];



memset(strCRegion, '\0' , strlen(seq)+1 ) ;



strcpy(strCRegion,seq);



}





int PowerSpectrum::doPowerSpectrumAnalysis() 

{





const double PI = 3.1415926535897932384626433832795;

//public static final double PI ;





complex<double>  i1(0.0,1.0);

//complex<double>  i1 = sqrt(-1);



long lenStr = strlen(strCRegion);



delete []s;

delete []x;





cout<<"Coding Region Length:"<<lenStr<<endl;





s	= new double[lenStr/2 + 1];

x	= new double[lenStr/2 + 1];



//temporary check; 

x1  =  new double[lenStr/2 + 1];

xp  =  new double[lenStr/2 + 1];

xp1  =  new double[lenStr/2 + 1];

xr  =  new double[lenStr/2 + 1];

xr1  =  new double[lenStr/2 + 1];

 

double *sf  =  new  double[lenStr/2 + 1]; 



char nuclt[5]= "ATGC";

int const lenWindow = (lenStr/3)*3;

//int const lenWindow = lenStr;





double winsq		= (1.0/lenWindow)+(1.0/(lenWindow*lenWindow));

long  wincube	    = lenWindow*lenWindow*lenWindow;





int countA =0;

int countC =0;

int countG =0;

int countT= 0;



complex<double> s1(0.0,0.0);

double ss2=0.0;

double sq=0.0;

double s2=0.0;

 //temporary check

double p1=0.0;

double p2=0.0;

double p3=0.0;



double sbar =0.0;

double p=0.0;



	

//ofstream out("power");



	for(long i =0; i < lenWindow; i++ )

	{



	if( strCRegion[i] == 'A') 

		{ countA++; continue;}

  	if( strCRegion[i] == 'G') 

		{ countG++; continue;}

  	if( strCRegion[i] == 'C') 

		{ countC++; continue;}

  	if( strCRegion[i] == 'T') 

		{ countT++; continue;}



	} 



//out<<"Count A  :"<<countA <<"    countT :"<<countT <<"countG   :"<<countG <<"countC   :"<<countC <<endl;



 

sq = 0.0+ pow(countA,2)+pow(countG,2) + pow(countT,2) + pow(countC,2);



double f =0.0;



for(long k = 1 ; k <=lenWindow/2  ; k++ )

{



	s2		=0.0;

	sbar	=0.0;

	p		=0.0;

	f= ((1.0)*k)/lenWindow;



// for each Nucleotide. A T G C





	for(int i = 0; i <4 ; i ++ )

	{

		s1= complex<double>(0.0,0.0);





		for (int bi=0; bi<lenWindow; bi++)

		{

				if(strCRegion[bi] == nuclt[i] )

					s1= s1+exp(i1*2.0*PI*f*(bi+1.0));

		}



		s2 = s2+ (abs(s1)*abs(s1));

		//ss2=ss2  + (1/lenWindow*lenWindow)*(abs(s1)*abs(s1));



	}



	//store the value of sf for each freq.

	sf[k-1]=s2;



//	to calculater power specturm and code region 



		sq   = sq/wincube;

		sbar =(winsq-sq);

		sq   = sbar * (lenWindow/2.0);



// to calcuater signal to noise ratio



		s2 = s2/(sq*lenWindow*lenWindow);



		p= s2/sbar;

  //temporary check;

        //p1=(p*lenWindow*sbar)/2;



	

		//p2=p1*2;

		p2=p/(lenWindow*sbar);

		p3=2*p2;



	s[k-1]=f;

	x[k-1]=p;

	xp[k-1]=p2;

	xp1[k-1]=p3;





	//temporary check

	//x1[k-1]=p1;

	//x2[k-1]=p2;

	//x3[k-1]=p3;





//	out<<"\n Frequency:"<<f <<" Power  : "<<p;



} // outer for loop



//calculating the sbar value



 double sbar1=0.00;

 double sbar2=0.00;



  for(k = 1 ; k <=lenWindow/2  ; k++ )

  {

	  sbar1=sbar1+sf[k-1];



  }

  sbar1=(sbar1*2)/lenWindow;

  sbar2=sbar1/(lenWindow*lenWindow);

  //calculating the sf/sbar1

  



  for(k = 1 ; k <=lenWindow/2  ; k++ )

  {

	   x1[k-1]=sf[k-1]/sbar1;



	  xr[k-1]=x1[k-1]/(lenWindow*sbar2);

	  

	  xr1[k-1]=xr[k-1]*2;



  }







int numPeaks=0;

int genePeak=0;



for( k = 1 ; k <= (lenWindow/2)  ; k++ )

{



	if( x[k-1] > 4.0)

	{

	   //if(  s[k-1] == 1.0/3.0  )] //added by Sreenu on mandars suggestion on 19/04/2002

	 if( ( s[k-1] >= 0.32 ) && ( s[k-1] <= 0.34 ))

             genePeak=1;

	   else

		   numPeaks++;

	}



}





//delete []s	  ;

//delete []x	  ;

//s =NULL  ;

//x =NULL  ;





if(genePeak==1)

{  

	if(numPeaks==0)

		 return 1;

	 else

		 return 2;

}

else

{

	return 0;

}







}  // end of member function.











//End of the pwerspectrum module







































//Included the Windowanalysis module



struct structCR

{

 long start;

 long end;

 char type;

};





class WindowAnalysis //__declspec(dllexport)   

{

public:



		

	long    *s;

	double   *x; 



	vector<structCR> vecCRegion;

			 WindowAnalysis();

			 WindowAnalysis(int len ,double thold,long intron , long exon);

	virtual  ~WindowAnalysis();

	void	  setWindowLen(int len);

	void	  settherHold(double thold);

	void      setIntronExon(long intron , long exon );

	int		  doFt(const char  * strSeq , char *logfile);

    long    framecount;



private :



	int		lenWindow;

	double  thersHold;

	long	exonCutOff;

	long	intronCutOff;

	int		findexp(int WindowSize);

	//long    framecount;

	complex<double>  * expval;

	WindowAnalysis& WindowAnalysis::operator =(const WindowAnalysis &cwa);

	WindowAnalysis(const WindowAnalysis &cwa);



};



















WindowAnalysis::WindowAnalysis()

{

	expval = NULL;

	 s	   = NULL;

     x	   = NULL;

	 framecount=0;





}



WindowAnalysis::~WindowAnalysis()

{

	

	delete []expval;	

	delete []s	  ;

	delete []x	  ;

	printf("Destructor of WindowAnalysis ");

	fflush(stdout);

  //delete []CRegion;





}



WindowAnalysis::WindowAnalysis(int len , double thold ,long intron , long exon)

{



	lenWindow=len;

	thersHold=thold;

	exonCutOff=exon;

	intronCutOff =intron;

	

	s	   = NULL;

    x	   = NULL;



	framecount=0;

//	cout<<"ThresHold Value "<<thersHold;

//	cout<<"Intron :"<<intronCutOff<<" Exon :"<<exonCutOff<<endl;



    expval=new complex<double>[lenWindow];



}





void WindowAnalysis::settherHold(double thold)

{

	thersHold =thold;

//	 cout<< "ThresHold Value "<<thersHold;



}



void WindowAnalysis::setWindowLen(int len)

{

    lenWindow=len;

	delete [] expval;

	expval=new complex<double>[lenWindow];



}

void  WindowAnalysis::setIntronExon(long intron , long exon )

{

	 exonCutOff=exon;

	 intronCutOff =intron;





}

int WindowAnalysis::doFt(const char  * strSeq , char * logfile)

{





const double F3 = 1.000/3;

const double PI = 3.14159265;





//ofstream	out("output" );

ofstream	out(logfile ,ios::app);



long lenStr = strlen(strSeq);



if(lenStr<lenWindow)

{

	out<<"Sequence length is less than windowlength."<<endl;

	return 1;

}



cout<<"\n Sequence Length :"<<lenStr<<endl;





int result= findexp(lenWindow);



//long *CRegion = new long[lenStr-lenWindow+ 1];

//long *s	    = new long[lenStr-lenWindow+ 1];

//double *x	= new double[lenStr-lenWindow+1];



delete [] s;

delete [] x;

s	    = new long[lenStr-lenWindow+ 1];

x	= new double[lenStr-lenWindow+1];





char nuclt[5]= "ATGC";



double winsq		= (1.0/lenWindow)+(1.0/(lenWindow*lenWindow));

long  wincube	    = lenWindow*lenWindow*lenWindow;





int countA =0;

int countC =0;

int countG =0;

int countT= 0;



complex<double> s1(0.0,0.0);



double sq=0.0;

double s2=0.0;

double sbar =0.0;

double p=0.0;



framecount=0;







for(long k = 0 ; k < lenStr-lenWindow+1 ; k++ )



{

	sq		=0.0;

	s2		=0.0;

	sbar	=0.0;

	p		=0.0;



	// Frame wise  count of A C G T seperately





if(k==0)

{

	for(long i =k; i < k+lenWindow; i++ )

	{



	if( strSeq[i] == 'A') 

		{ countA++; continue;}

  	if( strSeq[i] == 'G') 

		{ countG++; continue;}

  	if( strSeq[i] == 'C') 

		{ countC++; continue;}

  	if( strSeq[i] == 'T') 

		{ countT++; continue;}



	} //end of count for look.



}else

{



switch(strSeq[k-1])

{

case 'A' :

    { countA--; break;	}

case 'T' :

    { countT--; break;	}

case 'G' :

    { countG--; break;	}

case 'C' :

    { countC--; break;	}

}



switch(strSeq[k+lenWindow-1])

{

case 'A' :

    { countA++; break;	}

case 'T' :

    { countT++; break;	}

case 'G' :

    { countG++; break;	}

case 'C' :

    { countC++; break;	}

}



 

}



//out<<"Count A  :"<<countA <<" Frame count :"<<k<<endl;

//out<<"Count T  :"<<countT <<" Frame count :"<<k<<endl;

//out<<"Count G  :"<<countG <<" Frame count :"<<k<<endl;

//out<<"Count C  :"<<countC <<" Frame count :"<<k<<endl;



 sq = 0.0+ pow(countA,2)+pow(countG,2) + pow(countT,2) + pow(countC,2);



// for each Nucleotide. A T G C





	for(int i = 0; i <4 ; i ++ )

	{

		s1= complex<double>(0.0,0.0);





		for (int bi=0; bi<lenWindow; bi++)

		{

				if(strSeq[k+bi] == nuclt[i] )

					s1= s1+expval[bi];

		}



		s2 = s2+ (abs(s1)*abs(s1));



	}





//	to calculater power peak 



		sq   = sq/wincube;

		sbar =(winsq-sq);

		sq   = sbar * (lenWindow/2.0);



// to calcuater signal to noise ratio



		s2 = s2/(sq*lenWindow*lenWindow);



		p= s2/sbar;





	s[framecount]=k+(lenWindow/2);

	x[framecount]=p;



	out<<s[framecount]<<"                "<<x[framecount]<<"\n";

    framecount++;





} // End of Outer for loop(Frame wise)



//return(0);

//} 

// end of Ft function..

//void WindowAnalysis::findCodingRegion(long seqfileId , const char *  seqId)

//{



long *CRegion = new long[framecount];

//memset(CRegion , 0 ,(framecount+1)*sizeof(long));



int    flag=0; 

long i =0 , m =0;



	for (m =0; m<framecount; m++)

	{

       

		if(x[m]>=thersHold)

		{

			if(flag == 0)

			{

				

				CRegion[i] =s[m];

				//out<< "\n START :" <<s[m]<<" : "<<x[m];

				//i++;

				flag =1;

			}



		}



        else

        {

            if(flag==1)

       	     {

               if(x[m]<thersHold)

               {

				 	

					 	 out<< "\n START :" <<CRegion[i];

						 out<<"\t\tEND: "<<s[m-1]<<" : "<<x[m-1];		 

						i++;

						CRegion[i]=s[m-1];

						i++;

				}

				flag=0;



			 }

			

        } // end of if.

	}// End Coding region loop. 







if (flag==1)

{

	out<< "\n START :" <<CRegion[i];

	out<<"\t\tEND: "<<s[m-1]<<" : "<<x[m-1];		 

	i++;

	CRegion[i]=s[m-1];

	i++;

}



//delete []CRegion;

// that data is deleted 07/02/02 To inprove performance.



//delete []s	  ;

//delete []x	  ;

//s=NULL;

//x=NULL;

 

out<<"\n Final Coding Regions   :"<<endl;



flag =0;

int code =0;

long sta=0 , stp =0;

structCR obj;

obj.type=' ';

					

for(int j=0 ; j<i ; j +=2)

   {

	if(flag ==0)

	{

		sta = CRegion[j];

		stp = CRegion[j+1];

        if ( (stp-sta)>= exonCutOff ) code =1; else code =0; 	

		flag = 1 ;

    }

	else

	{                      

		if( (CRegion[j]-CRegion[j-1]) >=intronCutOff  )

		{

			if(code ==1)

			{	

				out <<"\nStart :"<<sta <<"  End  :" << stp;  			

				// inset in database  start new set.

				    obj.start=sta;

					obj.end=stp;

					vecCRegion.push_back(obj);

				//insetcodes( sta , stp,seqfileId,seqId);

				sta = CRegion[j];

				stp = CRegion[j+1];

			

				if ( (stp-sta)>= exonCutOff ) code =1; else code =0; 	

				flag = 1 ;

			}

			else

			{ 

				//cout << forget the old set.

				  if ((stp-sta)>= exonCutOff ) 

				  { out <<"\nStart :"<<sta <<"  End  :" << stp;

					obj.start=sta;

					obj.end=stp;

					vecCRegion.push_back(obj);

					//insetcodes( sta , stp,seqfileId,seqId);

				  }

				// may need to be  commented later.

				sta = CRegion[j];

				stp = CRegion[j+1];

				if ( (stp-sta)>= exonCutOff ) code =1; else code =0; 	

				flag = 1 ;

			}

		}

		else

		{



        	if( ( CRegion[j+1]-CRegion[j] ) >= exonCutOff ) code =1;

			stp = CRegion[j+1];

        	flag =1;



		}

	}

 }



if ((code == 1) || ((stp-sta) >exonCutOff) )

{



	out <<"\nStart :"<<sta <<"  End  :" << stp;

	obj.start=sta;

	obj.end=stp;

	vecCRegion.push_back(obj);

	//mapCRegion[sta] =stp;

 	//insetcodes( sta , stp,seqfileId,seqId);		

}





delete []CRegion; 

CRegion=NULL;





vector<structCR>::iterator it;



PowerSpectrum powerObj;



char * strtemp=NULL;



for(it =vecCRegion.begin(); it !=vecCRegion.end(); it++)

{

	//long start =it->start;

	//long end =it->end;



	strtemp = new char[( (it->end - it->start) +1)];



	memset(strtemp , '\0' ,(it->end - it->start) +1 );



	strncpy( strtemp , &strSeq[it->start] ,(it->end - it->start) );



	//strncpy( strtemp , &strSeq[start] , end -start);



    printf("Length : %d \n ", strlen(strtemp));



	printf("start %ld , End  %d  length %d\n " , it->start ,it->end , (it->end - it->start));

	fflush(stdout);



	powerObj.setSequence(strtemp);

	int peaks =powerObj.doPowerSpectrumAnalysis();

	if(peaks==0)

	{

		it=vecCRegion.erase(it);

		if( it !=vecCRegion.end() )

		--it;

		else

			break;

	}

	else if(peaks==1)

	{

		it->type ='S';

	}

	else if(peaks==2)

	{

		it->type ='M';

	}

delete []strtemp;

strtemp = NULL;

}







  //store the coding regions in an array





//int elsize;

 //elsize=vecCRegion.size();







delete []strtemp;

strtemp = NULL;





//delete []CRegion; 

//CRegion=NULL;



//delete []s	  ;

//delete []x	  ;



return(0); 



} // end of Coding region finding.

 



// to find exponential positional values to window size .



int WindowAnalysis::findexp(int WindowSize)

{





complex<double>  i1(0.0,1.0);

const double PI = 3.14159274;

const double F3 = 1.000/3;



for(int k=0; k<WindowSize; k++)

{ 

	

	expval[k] = exp(-i1*2.0*PI*F3*(k+1.0));		



}



return(0);



}



// ///////////////// End //////////////







//End of Windowanalysis























struct storeseq

{

	char seqhead[100];



	storeseq(char *s)

	{strcpy(seqhead,s);}



};



struct storeth

{

	char seqt[100];

	double sf;

	storeth(char *s,double b)

	{strcpy(seqt,s);sf=b;}





};



 class Chpcity

 {



	 

   

        



   FILE *fp1,*fp2,*fp3,*fp4;

   char *seq;

   fpos_t begpos,endpos;

   vector <storeseq> stseq;

   vector <storeseq>::iterator it;

   

   vector <storeth> stth;

   vector <storeth> stth1;

   vector <storeth> ::iterator it2;

   vector <storeth> ::iterator it3;









 public:

	 void doGenscan();

	 void doBlast();

	 void doParse();

     void dofft();

	 void filterSeq();





      ~Chpcity()

	  {

		 cout<<"Destructor of chpcity"<<endl;

		 stseq.clear();

	  }





 };

 



 void Chpcity::doGenscan()

 {

	 

	 //variable where the sequence will be stored



	 WindowAnalysis winobj;

     int seqno=1;

	 //storeth objth;



	 

	 char fname[60];

	 char ch,ch1;

	 int sin;

	 char seqhead[100];



	 vector<structCR>::iterator it1;



	 

	 winobj.settherHold(0.0);

	 winobj.setIntronExon(50,50);

	 

   

	 



	 printf("Enter the file name:");

	 gets(fname);

	 

	 fp1=fopen(fname,"r");



	 	 

	 if(fp1==NULL)

	 {

		 cout<<"File doesn't exist"<<endl;

		 exit(1);

	 }



	 //create a temporary file to store the sequence without the newline characters

	 

	 //cout<<"Enter the file name for non coding sequences:");

	 //gets(fname);

	 

	 fp3=fopen("d:\\tholdamy.txt","w");

	 

	 ch=fgetc(fp1);

	 

	 while(ch!=EOF)

	 {

		 sin=0;

		

	    fp2=fopen("tmpseq.tmp","w+");

	 

	    if(fp2==NULL)

		{

		    cout<<"Could not create temporary file"<<endl;

		    exit(1);

		}

			 

		do 

		{

		 

		 //check for the initial '>'symbol put the characters after the next line

		 if(ch=='>')

		 {

		   while(true)

		   {

		    ch=fgetc(fp1);

			 

			 if(ch=='\n')

			 {

				 ch=fgetc(fp1);

				 break;

			 }

			 seqhead[sin]=ch;

			 sin++;

		   }//end ofwhile	

		   seqhead[sin]='\0';

		 



		  }//end of if

		



		 if(ch!='\n')  

		   fputc(ch,fp2);



		 ch=fgetc(fp1);



		 if(ch==EOF)

		 {

		   break;

		 }

	 		

		}while(ch!='>');//end of while





		



        fgetpos(fp2,&endpos);



		//int endpos=begpos;

		

		seq=new char[endpos +2];

		memset(seq,'\0',endpos +2);





		begpos=0;

		fsetpos(fp2,&begpos);



		ch1=fgetc(fp2);

		int seqin=0;



		while(ch1!=EOF)

		{

			seq[seqin]=toupper(ch1);

			seqin++;

			ch1=fgetc(fp2);

		}





		/*char exseqh[30];



		for(int j1=0;seqhead[j1]!=',';j1++);



		for(int ie=j1+1,i=0;seqhead[ie]!='\0';ie++,i++)

		{

			exseqh[i]=seqhead[ie];

		}



		exseqh[i]='\0';*/



		//if(strcmp(seqhead,"(gi|16127994:1080570-1080686, 1080677-1081408), b1016m")==0)

		//{

		//winobj.setWindowLen(strlen(seq)); 

			

		winobj.setWindowLen(300); 

		winobj.doFt(seq,"log");



		//cout<<winobj.x[0]<<endl;



		//objth.sf=winobj.x[0];

		

		//strcpy(objth.seqt,seqhead);

		//double tx=winobj.x[0];



		//cout<<objth.sf<<endl;

		

		for(int fr=0;fr<winobj.framecount;fr++)

		  stth.push_back(storeth(seqhead,winobj.x[fr]));

		//}

		 

		 	

	     fclose(fp2);

		

		remove("tmpseq.tmp");

        winobj.vecCRegion.clear();    

		seqno++;

		

		delete []seq;

      



	  }

	  char tempseqt[100];

	  float tempsf;

	  



	  for(it2=stth.begin();it2!=stth.end();it2++)

	  {		   		  

		  it3=it2;

		  it3++;

		  for(;it3!=stth.end();it3++)

		  {

			  if(it2->sf < it3->sf)

			  {

				  tempsf=it3->sf;

				  it3->sf=it2->sf;

				  it2->sf=tempsf;



				  strcpy(tempseqt,it3->seqt);

				  strcpy(it3->seqt,it2->seqt);

				  strcpy(it2->seqt,tempseqt);

				  



			  }

		  }

		  

	  }



	   

	  fprintf(fp3,"%s\n",fname); 

	  for(it2=stth.begin();it2!=stth.end();it2++)

	  {

		  fprintf(fp3,"%0.2f\n",it2->sf);

		  //cout<<it2->seqt<<" "<<it2->sf<<endl;

	  }



       



	  

	  _fcloseall();



	  //system("del d:\\tmpseq.tmp");





 

 

 

 

 }



 

 

 void Chpcity::dofft()

 {

	 

	 PowerSpectrum powobj;



	 char heading[150];

	 int hin;

	 char filename[50];  

	 char ch,ch1;

	 int lenseq;



	 int no=1;

	 int no1=1;





	 printf("Enter the file name for powerspectrum anal.:");

	 gets(filename);



	 fp1=fopen(filename,"r");

	 

	 if(fp1==NULL)

	 {

		 cout<<"Cannot Open file"<<endl;

		 exit(1);

	 }



	 

     



	 //fp4=fopen("d:\\anno\\ecolsfr.txt","w");



    

	 ch=fgetc(fp1);

	 

	 while(ch!=EOF)

	 {

		 hin=0;

		

	    fp2=fopen("tmpseq.tmp","w+");

	 

	    if(fp2==NULL)

		{

		    cout<<"Could not create temporary file"<<endl;

		    exit(1);

		}

			 

		do 

		{

		 

		 //check for the initial '>'symbol put the characters after the next line

		 if(ch=='>')

		 {

		   while(true)

		   {

		    ch=fgetc(fp1);

			 

			 if(ch=='\n')

			 {

				 ch=fgetc(fp1);

				 break;

			 }

			 

			 heading[hin]=ch;

			 hin++;

		   }//end ofwhile	



		   heading[hin]='\0';

		 

		  }//end of if

		

		 



		  if(ch!='\n')  

		   fputc(ch,fp2);



		  ch=fgetc(fp1);



		 if(ch==EOF)

		 {

		   break;

		 }

	 		

		}while(ch!='>');//end of while





        fgetpos(fp2,&endpos);



		//int endpos=begpos;

		

		seq=new char[endpos +2];

		memset(seq,'\0',endpos +2);





		begpos=0;

		fsetpos(fp2,&begpos);



		ch1=fgetc(fp2);

		int seqin=0;



		while(ch1!=EOF)

		{

			seq[seqin]=toupper(ch1);

			seqin++;

			ch1=fgetc(fp2);

		}



		 //cout<<seq<<endl;

		char exseqh[50];



		for(int j1=0;heading[j1]!=',';j1++);



		for(int ie=j1+1,i=0;heading[ie]!='\0';ie++,i++)

		{

			exseqh[i]=heading[ie];

		}



		exseqh[i]='\0';

		

				

		

		

		fp3=fopen(exseqh,"w");



		if(fp3==NULL)

		{

			cout<<"Could not create File"<<endl;

		}



		powobj.setSequence(seq);

		int peak=powobj.doPowerSpectrumAnalysis();



		lenseq=strlen(powobj.strCRegion);

		

		lenseq=(lenseq/3)*3;



		//fprintf(fp3,">%s\n",heading);



		//fprintf(fp4,">%s\n",exseqh);



        //fprintf(fp3,"%s\n",seq);



        

		char cfreq[7];

		float ffreq;



	 	 //double frq;

	  	 //frq=1.0/3.0;

		

		for(long k1=0;k1<lenseq/2;k1++)

		{

			



			//if(powobj.s[k1]==frq)

			//{



				//sprintf(cfreq,"%0.4f",powobj.s[k1]);

				//ffreq=atof(cfreq);

				//ffreq=1/ffreq;

				//cout<<powobj.s[k1]<<" "<<powobj.x1[k1]<<" "<<powobj.x[k1]<<endl;				

				//if(powobj.x1[k1]>4.0)

				//{

			       //stth.push_back(storeth(heading,powobj.x1[k1]));

				//}

				/*else

				{

					

		            stth1.push_back(storeth(exseqh,powobj.x1[k1]));

				}*/



	

			  fprintf(fp3,"%f %f\n",powobj.s[k1],powobj.x[k1]); 



			  no++;

			//}



		}

		







		//stth.clear();

		//stth1.clear();



		//fprintf(fp3,"\n");

         

		fclose(fp2);



		remove("tmpseq.tmp");

		delete []seq;





   

		

     }//end of while



	   /*char tempseqt[100];

		double tempsf;

		

		for(it2=stth.begin();it2!=stth.end();it2++)

		{		   		  

		  it3=it2;

		  it3++;

		  for(;it3!=stth.end();it3++)

		  {

			  if(it2->sf < it3->sf)

			  {

				  tempsf=it3->sf;

				  it3->sf=it2->sf;

				  it2->sf=tempsf;



				  strcpy(tempseqt,it3->seqt);

				  strcpy(it3->seqt,it2->seqt);

				  strcpy(it2->seqt,tempseqt);

				  



			  }

		  }

		  

		}



		

		



		for(it2=stth.begin();it2!=stth.end();it2++)

		{

		  fprintf(fp3,"%d %s %0.14lf\n",no,it2->seqt,it2->sf);

		  no++;

		}*/









		/*for(it2=stth1.begin();it2!=stth1.end();it2++)

		{		   		  

		  it3=it2;

		  it3++;

		  for(;it3!=stth1.end();it3++)

		  {

			  if(it2->sf < it3->sf)

			  {

				  tempsf=it3->sf;

				  it3->sf=it2->sf;

				  it2->sf=tempsf;



				  strcpy(tempseqt,it3->seqt);

				  strcpy(it3->seqt,it2->seqt);

				  strcpy(it2->seqt,tempseqt);

				  



			  }

		  }

		  

		}



		

		for(it2=stth1.begin();it2!=stth1.end();it2++)

		{

		  fprintf(fp4,"%d %s %0.14lf\n",no1,it2->seqt,it2->sf);

		  no1++;

		}*/





 _fcloseall();





 }



 

 void Chpcity::doParse()

 {



     char pfname[50];

	 char chp;



	 //printf("Enter the file name to parse:");

	 //gets(pfname);



	 strcpy(pfname,"c:\\localb\\testbl.txt");



	 fp1=fopen(pfname,"r");



	 if(fp1==NULL)

	 {

		 cout<<"File dooesn't exist"<<endl;

		 exit(1);

	 }



	 chp=fgetc(fp1);



	 char *seqsrch,*posit;

	 

	 



	 char seqname[30];

	 

 	 while(chp!=EOF)

	 {



	   while(chp!=EOF)

	   {

		   seqsrch=new char[7];

	       memset(seqsrch,'\0',7);

		  



		  for(int i=0;i<6;i++)

		  {

			 		 

			if(chp=='\n' || chp==' ')

			{

				chp=fgetc(fp1);

				break;

			}



			 seqsrch[i]=chp;

			 chp=fgetc(fp1);

			 

			 if(chp==EOF)

				 break;



			

		  }

		   //seqsrch[i]='\0';



		   if(strcmp(seqsrch,"Query=")==0)

			  break;

		   

		   delete []seqsrch;

	   }

	        

        

	         if(chp==EOF)

				 break;

		



	 

	 //extract the sequence name

		 int j=0;

		 chp=fgetc(fp1);

		



		 while(chp!='\n')

		 {

			 

			  seqname[j]=chp;

			  j++;

			  chp=fgetc(fp1);





		 }

		 seqname[j]='\0';



	 //check for hits and extract positives



		 chp=fgetc(fp1);



		int nohit=0;

		int nomatch=0;

		char positptage[4];//store the percentage from the positives



		while(chp!=EOF)

		{

		   posit=new char[10];

	       memset(posit,'\0',10);



		   for(int i=0;i<9;i++)

		   {

			 						   

			 if(chp=='\n' || chp==' ')

			 {

				 chp=fgetc(fp1); 

				 break;

			 }



			 

			 if(chp=='*')//*means there is no hit for the sequence

			 {

				 nohit=1;

				 break;

			 }

			

			 posit[i]=chp;

			 chp=fgetc(fp1);



			 if(strcmp(posit,"Lambda")==0)

			 {

				 nomatch=1;

				 break;



			 }

						 

						

		   }

		    

		    //posit[i]='\0';



			if(nohit==1)//no hits found for the sequence

				break;

			

			if(nomatch==1)//no positives >=60%

				break;

			

		    if(strcmp(posit,"Positives")==0 )

			{

			   //check for >=60%

				

				while(chp!='(')

				{

					chp=fgetc(fp1);

				}



				chp=fgetc(fp1);

				

				int pspin=0;

				

				while(chp!='%')

				{

					positptage[pspin]=chp;

					pspin++;

					chp=fgetc(fp1);

				}

				

				positptage[pspin]='\0';



				int intptage=atoi(positptage);



				if(intptage>=60)

				{

                    stseq.push_back(storeseq(seqname));

					break;

				}

			

			}//end of if

			delete []posit;



		}

		

	}

   

	fclose(fp1);

	

	

	for(it=stseq.begin();it!=stseq.end();it++)

	{

		cout<<it->seqhead<<endl;



	}



 

	

 }



 void Chpcity::filterSeq()

 {

	 //this module filters out the sequences that has Positives <=x%



	 char fch;

	 char seqhead[30];

	 int findseq,fseqin;



	 fp1=fopen("d:\\noncodseq\\noncoseq.txt","r");

	 

	 if(fp1==NULL)

	 {

		 cout<<"File doesn't exist"<<endl;

		 exit(1);



	 }





	 fp2=fopen("d:\\noncodseq\\fnoncoseq.txt","w");



	 

	 fch=fgetc(fp1);



	for(it=stseq.begin();it!=stseq.end();it++)

	{

		findseq=0;		 

	 

	  while(findseq==0)

	  {

		  fseqin=0;

	

		 //check for the initial '>'symbol put the characters after the next line

		 if(fch=='>')

		 {

		   while(true)

		   {

		    fch=fgetc(fp1);

			 

			 if(fch=='\n')

			 {

				 fch=fgetc(fp1);

				 break;

			 }

			 

			 seqhead[fseqin]=fch;

			 fseqin++;

		   }//end ofwhile	



		   seqhead[fseqin]='\0';

		 

		  }//end of if



		 if(strcmp(it->seqhead,seqhead)==0)

		 {

			 fprintf(fp2,">%s\n",seqhead);

			 while(fch!='\n')

			 {

				 fputc(fch,fp2);

				 fch=fgetc(fp1);



			 }

			 fputc('\n',fp2);



			 findseq=1;



		 }



          fch=fgetc(fp1);

		



	  }

	 

    

	}//end of for loop



  

	_fcloseall();

 

 }



 

 void Chpcity::doBlast()

 {

	 	 

	 //execute the blastall program

	 int retval=system("c:\\localb\\blastall -p blastx -d c:\\localb\\pdbaa -i d:\\noncodseq\\noncoseq.txt -o c:\\localb\\testbl.txt -e 0.1");

	  

 }



 

 

 int main(int argc, char* argv[])

 {

	Chpcity Ch;

	//Ch.doGenscan();

	//Ch.doBlast();

	//Ch.doParse();

	

	//Ch.filterSeq();

	Ch.dofft();

	

    //{//blastall -p blastx -d pdbaa -i d:\noncoseq.txt -o test1.txt -e 0.1}

    //blastall -p blastx -d pdbaa -i d:\noncoseq.txt -o test.xml -m 7 -e 0.1



	

	return 0;

 }