//file:combine.cc
//Qing Chang
//qchan@sac.uky.edu
//last modified on 11-20-98
//project: prediction of transmembrane protein

/*************************************************************************
input:amino acid string, the single letter version
output:the sequence predicted as helix or sheet, in the single letter
*************************************************************************/

#include <iostream.h>
#include <string.h>


#define TRUE 1
#define FALSE 0
#define INPEP 650       //the maximum length of the input polypeptide
#define H_CORE 6        //the length of the core region of helix
#define S_CORE 5        //the length of the core region of sheet    
#define EXTENSION 4     //the length of the extended region
#define key_helix 1.15  //the boarder line of helix prediction
#define key_sheet 1.2   //the boarder line of sheet prediction

#define LH 1.373331  /*************P-helix value for each aa***********/
#define IH 1.264827
#define FH 1.201535
#define VH 1.194153
#define AH 1.179271
#define MH 1.13655
#define WH 1.115379
#define RH 1.036602
#define YH 1.017381
#define CH 1.01368
#define TH 0.965942
#define SH 0.928515
#define HH 0.906637
#define KH 0.894321
#define EH 0.806584
#define GH 0.772239
#define QH 0.740862
#define DH 0.728439
#define NH 0.694558
#define PH 0.418516


#define CS 1.811841 /************P-sheet value for each aa************/
#define YS 1.363843
#define DS 1.362563
#define QS 1.269034
#define WS 1.226841
#define NS 1.221741
#define VS 1.219667
#define GS 1.159773
#define TS 1.159232
#define RS 1.048056
#define KS 0.955804
#define SS 0.940886
#define MS 0.918993
#define IS 0.91475
#define ES 0.871015
#define FS 0.815549
#define LS 0.769492
#define PS 0.732135
#define HS 0.682322
#define AS 0.665627

void instruction ();//display instruction on the screen
char* get_input ();//get input
void prediction (char*, int);//translate peptide to double array
void helix_core_prediction (double*, int,int&, int&);//predict helix core
void helix_extension(double*, int&, int&, int);//predict helix extension 
void helix_prediction (double*, int, char*, char*);//predict helix structure 
void sheet_prediction (double*, int, char*, char*);//predict sheet structure
void sheet_core_prediction (double*, int,int&, int&);//predict sheet core
void sheet_extension(double*, int&, int&, int);//predict sheet extension


int main()
{
   char* inpep; //user input polypeptide string

   instruction ();//display the instruction

   inpep = new char[INPEP];
   inpep = get_input ();//get user's input
   prediction (inpep, INPEP);//predict structure 
   delete inpep;
   return 0;
}

void instruction ()
{
   cout <<endl
	<<"*********************INSTRUCTION*****************************"<<endl
	<<"This program is to predict the 2nd structure of transmembrane"<<endl
	<<"proteins. Please type in or paste your polypeptide, one letter"<<endl 
	<<"per amino acid. Then hit the return. The predicted helix or" <<endl
	<<"sheet region will be displayed as results." <<endl
	<<"*************************************************************"<<endl
	<<endl;
   return;
}


char* get_input ()
{
   char* local;//local string pointer
   local = new char[INPEP];

   cout <<"Please type in or paste the polypeptide you like to predict:";
   cin.getline (local, INPEP,'*');

   return local;
}

void prediction (char* peptide, int size)
{
    double p_value[size];//array holds the p_value of each aa in the peptide array
    double s_value[size];//array holds the sheet value of each aa
    char result[size]; //array holds prediction result
    for (int i=0; i<size; i++)
       result[i]='.';

    int real_size=0;//the real size of the input peptide
    int j=0;
    for (int i=0; i<size; i++)      /* assign value to p_value array */ 
    {
       switch ( peptide[i] )        
       {
	  case 'L':
	  case 'l':
	     p_value[j] = LH;
	     s_value[j] = LS; j++; break;
          case 'I':
	  case 'i':
	     p_value[j] = IH;
	     s_value[j] = IS; j++; break;
          case 'F':
	  case 'f':
	     p_value[j] = FH;
	     s_value[j] = FS; j++; break;
          case 'V':
	  case 'v':
	     p_value[j] = VH;
	     s_value[j] = VS; j++; break;
          case 'A':
	  case 'a':
	     p_value[j] = AH;
	     s_value[j] = AS; j++; break;
          case 'M':
	  case 'm':
	     p_value[j] = MH;
	     s_value[j] = MS; j++; break;
          case 'W':
	  case 'w':
	     p_value[j] = WH;
	     s_value[j] = WS; j++; break;
          case 'R':
	  case 'r':
	     p_value[j] = RH;
	     s_value[j] = RS; j++; break;
          case 'Y':
	  case 'y':
       	     p_value[j] = YH;
	     s_value[j] = YS; j++; break;
          case 'C':
	  case 'c':
	     p_value[j] = CH;
	     s_value[j] = CS; j++; break;
          case 'T':
	  case 't':
	     p_value[j] = TH;
	     s_value[j] = TS; j++; break;
          case 'S':
	  case 's':
	     p_value[j] = SH;
	     s_value[j] = SS; j++; break;
          case 'H':
          case 'h':
	     p_value[j] = HH;
	     s_value[j] = HS; j++; break;
          case 'K':
	  case 'k':
	     p_value[j] = KH;
	     s_value[j] = KS; j++; break;
          case 'E':
	  case 'e':
	     p_value[j] = EH;
	     s_value[j] = ES; j++; break;
          case 'G':
	  case 'g':
	     p_value[j] = GH;
	     s_value[j] = GS; j++; break;
          case 'Q':
	  case 'q':
	     p_value[j] = QH;
	     s_value[j] = QS; j++; break;
          case 'D':
	  case 'd':
	     p_value[j] = DH;
	     s_value[j] = DS; j++; break;
          case 'N':
	  case 'n':
	     p_value[j] = NH;
	     s_value[j] = NS; j++; break;
          case 'P':
	  case 'p':
	     p_value[j] = PH;
	     s_value[j] = PS; j++; break;
          default:
       	    p_value[j] = -1;
	    s_value[j] = -1; j++; 
       }//end of switch
       if (p_value[j-1]==-1)
       {
	  real_size=j;
	  i=size;
       }
    }//end of for loop
    
    helix_prediction(p_value, real_size, peptide, result);
    sheet_prediction(s_value, real_size, peptide, result);
 

    for (int i=0; i<(real_size); i++)
    {
	if (result[i]=='H')
	   cout <<"The "<<i+1<<" residue is in Helix"<<endl;
        if (result[i]=='S')
	   cout<<"The "<<i+1<<" residue is in Sheet"<<endl;
    }

    return;
}


void helix_prediction (double* p_value, int real_size, char* peptide, char* result)
{
    int start =0;//start position of core prediction
    int end = start + H_CORE-1;//end position of core prediction
    int local_result = FALSE;//prediction result
    int extra=0;//the count of non amino acid in the input

    while (end<(real_size-1))
    {
       helix_core_prediction(p_value, start,end, local_result);
       if (local_result==TRUE) //predicted helix core
       {
	  helix_extension(p_value, start,end,real_size);
          int start2=start+1;
	  int end2=start2+H_CORE-1;
	  int result2=FALSE;
	  while ((start2<=end)&&(end2<(real_size-1))) //predict within the core
	  {
	     helix_core_prediction(p_value,start2,end2,result2);
	     if (result2==TRUE)
	     {
		helix_extension(p_value,start2,end2,real_size);
                if ((start2<start)&&(end2<=end))
		{
		   start=start2;
		   start2++; end2=start2+H_CORE-1; 
                }
		else if ((end2>end)&&(start2>=start))
		{
		   end=end2;
		   start2++;end2=start2+H_CORE-1;
                }
		else if ((start2<start)&&(end2>end))
		{
		   start=start2;end=end2;
		   start2++;end2=start2+H_CORE-1;
                }
		else
		{
		   start2++;
		   end2 = start2+H_CORE-1;
                }
            }//end of if 
            else
	    {
	       start2++;end2=start2+H_CORE-1;
            }
         }//end of while
	 for (int k=start;k<=end;k++)
            result[k]='H';
	 start = end+1; end=start+H_CORE-1;
      }//end of if         
       else
       {
	  start++;
	  end = start + H_CORE-1;
       }
    }//end of while

    return;
}

void sheet_prediction (double* s_value, int real_size, char* peptide, char* result)
{
   int start =0;//start position of core prediction
   int end = start + S_CORE-1;//end position of core prediction
   int local_result = FALSE;//prediction result
   int extra=0;//the count of non amino acid in the input

   while (end<(real_size-1))
     {
      sheet_core_prediction(s_value, start,end, local_result);
      if (local_result==TRUE) //predicted sheet core
      {
         sheet_extension(s_value, start,end,real_size);
	 int start2=start+1;
	 int end2=start2+S_CORE-1;
       	 int result2=FALSE;
	 while ((start2<=end)&&(end2<(real_size-1))) //predict within the core
         {
	     sheet_core_prediction(s_value,start2,end2,result2);
 	     if (result2==TRUE)
	     {
		sheet_extension(s_value,start2,end2,real_size);
		if ((start2<start)&&(end2<=end))
                {
	           start=start2;
		   start2++; end2=start2+S_CORE-1;
		}
		else if ((end2>end)&&(start2>=start))
		{
		   end=end2;
	           start2++;end2=start2+S_CORE-1;
		}
		else if ((start2<start)&&(end2>end))
	       {
		 start=start2;end=end2;
	         start2++;end2=start2+S_CORE-1;
	       }
               else
	       {
		  start2++;
	          end2 = start2+S_CORE-1;
	       }
	    }//end of if
	    else
	    {
	       start2++;end2=start2+S_CORE-1;
	    }
	 }//end of while
	 for (int k=start;k<=end;k++)
	    if (result[k]!='H')
	       result[k]='S';
	 start = end+1; end=start+S_CORE-1;
     }//end of if
     else
    {
       start++;
       end = start + S_CORE-1;
    }
  }//end of while
  return;
}

void helix_core_prediction (double* core, int start,int& end, int& result)
{
   result = FALSE;
   double sum = 0;
   for (int i=start; i<(start+H_CORE); i++)
   {
      sum = sum + core[i];
   }
   if ( (sum/H_CORE) >= key_helix )
   {
      result = TRUE;
   }

   return;
}

void sheet_core_prediction (double* core, int start,int& end, int& result)
{
   result = FALSE;
   double sum = 0;
   for (int i=start; i<(start+S_CORE); i++)
   {
      sum = sum + core[i];
   }
   if ( (sum/S_CORE) >= key_sheet )
  {
     result = TRUE;
  }

  return;
}


void helix_extension(double* extend, int& start, int& end, int size)
{
   double sum = 0;
   while(((start-EXTENSION)>=0)&&(sum!=-1)) /*entend to left*/
   {
      for (int i=start-1;i<=(start-EXTENSION);i--) 
      {
	 sum=sum+extend[i];
      }
      if (sum/EXTENSION >=key_helix)
         start = start- EXTENSION;
      else
         sum = -1;
   }

   sum =0;
   while(((end+EXTENSION)<(size-1))&&(sum!=-1))/*extend to right*/
   {
      for (int i=end+1;i<=(end+EXTENSION);i++)
      {
	sum = sum + extend[i];
      }
      if (sum/EXTENSION >= key_helix)
         end=end+EXTENSION;
      else 
	 sum=-1;
   }

   return;
}

void sheet_extension(double* extend, int& start, int& end, int size)
{
   double sum = 0;
   while(((start-EXTENSION)>=0)&&(sum!=-1)) /*entend to left*/
   {
       for (int i=start-1;i<=(start-EXTENSION);i--)
       {
          sum=sum+extend[i];
       }
       if (sum/EXTENSION >=key_sheet)
	  start = start- EXTENSION;
       else
	  sum = -1;
   }

   sum =0;
   while(((end+EXTENSION)<(size-1))&&(sum!=-1))/*extend to right*/
   {
      for (int i=end+1;i<=(end+EXTENSION);i++)
     {
         sum = sum + extend[i];
     }
     if (sum/EXTENSION >= key_sheet)
	 end=end+EXTENSION;
     else
	 sum=-1;
  }

  return;
}