/*  PAMxT.c       Jean-Michel Claverie      SGI     feb 1991 / March 1993 */
/* pamxt - Calculate a PAM to a given number;
*       - or/and a scoring matrix given a "transition" matrix
* -----------------------------------------------------------------------      
* The evolution model corresponding to this program is symetrical
*   [pam] [T] [pam]  with the default option,
*   asymetrical  [pam] [T]  with the -a (asymetrical) option
* ----------------------------------------------------------------------
* The  command syntax is
*
* 	pamxt <PAM distance>
*   or  
*       pamxt <PAM distance> <transition matrix> -a 
*       output file -> PAMxT.out
*
*---------------------------- IMPORTANT WARNING --------------------------
* The transition matrix must be given according the 
*          ARNDCQEGHILKMFPSTWYV   order !
*-------------------------------------------------------------------------
* Relevant literature:
*
* Claverie, J-M. (1993) J. MOL. BIOLOGY 234: 1140-1157.
* Detecting Frame Shifts by Amino Acid Sequence Comparison.
*
* Dayhoff, M. O., Schwartz, R. M. & Orcutt, B. C. (1978).
* A model of evolutionary change in proteins. In Atlas of Protein Sequence
* and Structure, Vol. 5, Suppl. 3, Ed. M. O. Dayhoff, pp. 345-352.
* Natl. Biomed. Res. Found., Washington.
*
* Karlin, S. & Altschul, S., Methods for assessing the statistical significance
* of molecular sequence features by using general scoring schemes
* Proc. Nat. Acad. Sci. USA, 87, 2264-2268, 1990
*/

/*  compile with  cc -lm     */
#include <stdio.h>
#include <math.h>
#define DIM 20
#define MAXLINE 1000
#define MLIMIT 20

/*	One letter amino acid code	*/

char aa[] = {"ARNDCQEGHILKMFPSTWYV"};

/*	Dayhoff amino acid background frequencies	*/

double fq[DIM] = {
0.086, 0.041, 0.040, 0.047, 0.033, 0.038, 0.050, 0.089, 0.034, 0.037,
0.085, 0.081, 0.015, 0.040, 0.051, 0.070, 0.058, 0.010, 0.030, 0.065};

/*	Dayhoff 1 PAM stochastic amino acid transition matrix	*/

int mdm[400] = {
9867,  2,  9, 10,  3,  8, 17, 21,  2,  6,  4,  2,  6,  2, 22, 35, 32, 0,  2, 18,
 1, 9914,  1,  0,  1, 10,  0,  0, 10,  3,  1, 19,  4,  1,  4,  6,  1, 8,  0,  1,
 4,  1, 9822, 36,  0,  4,  6,  6, 21,  3,  1, 13,  0,  1,  2, 20,  9, 1,  4,  1,
 6,  0, 42, 9859,  0,  6, 53,  6,  4,  1,  0,  3,  0,  0,  1,  5,  3, 0,  0,  1,
 1,  1,  0,  0, 9973,  0,  0,  0,  1,  1,  0,  0,  0,  0,  1,  5,  1, 0,  3,  2,
 3,  9,  4,  5,  0, 9876, 27,  1, 23,  1,  3,  6,  4,  0,  6,  2,  2, 0,  0,  1,
10,  0,  7, 56,  0, 35, 9865,  4,  2,  3,  1,  4,  1,  0,  3,  4,  2, 0,  1,  2,
21,  1, 12, 11,  1,  3,  7, 9935,  1,  0,  1,  2,  1,  1,  3, 21,  3, 0,  0,  5,
 1,  8, 18,  3,  1, 20,  1,  0, 9913,  0,  1,  1,  0,  2,  3,  1,  1, 1,  4,  1,
 2,  2,  3,  1,  2,  1,  2,  0,  0, 9871,  9,  2, 12,  7,  0,  1,  7, 0,  1, 33,
 3,  1,  3,  0,  0,  6,  1,  1,  4, 22, 9947,  2, 45, 13,  3,  1,  3, 4,  2, 15,
 2, 37, 25,  6,  0, 12,  7,  2,  2,  4,  1, 9924, 20,  0,  3,  8, 11, 0,  1,  1,
 1,  1,  0,  0,  0,  2,  0,  0,  0,  5,  8,  4, 9875,  1,  0,  1,  2, 0,  0,  4,
 1,  1,  1,  0,  0,  0,  0,  1,  2,  8,  6,  0,  4, 9944,  0,  2,  1, 3, 28,  0,
13,  5,  2,  1,  1,  8,  3,  2,  5,  1,  2,  2,  1,  1, 9924, 12,  4, 0,  0,  2,
28, 11, 34,  7, 11,  4,  6, 16,  2,  2,  1,  7,  4,  3, 17, 9840, 38, 5,  2,  2,
22,  2, 13,  4,  1,  3,  2,  2,  1, 11,  2,  8,  6,  1,  5, 32, 9869, 0,  2,  9,
 0,  2,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  1,  0,  1,  0, 9976, 1,  0,
 1,  0,  3,  0,  3,  0,  1,  0,  4,  1,  1,  0,  0, 21,  0,  1,  1, 2, 9947,  1,
13,  2,  1,  1,  3,  2,  2,  3,  3, 57, 11,  1, 17,  1,  3,  2, 10, 0, 2, 9901};


double	pam0[DIM][DIM];                             /* initial pam */
double  store[DIM][DIM];			    /* storage */
double  pamn [DIM][DIM];                            /* n th power of pam */
double T[DIM][DIM];                                 /* transition matrix */
double PM[DIM][DIM];                                /* freq*pamn matrix */
double PP[DIM][DIM];				    /* fq[i]*fq[j]      */
double MAT[DIM+3][DIM+3];		            /* Final Scoring matrix */

main(argc,argv)
int argc;
char *argv[];
{
int i,j,k,n,N,size, option=0;
double S, lambda;
FILE *F;

if (argc <2 || argc >4)
   { printf("\n USAGE   %s  number\n",argv[0]); 
     printf  ("         %s  number  Transition_Matrix  [-a]\n",argv[0]); exit(0);}

if (1 != sscanf(argv[1], "%d", &N))
   { printf("\n ERROR:  pam  number  %s\n",argv[1]); exit(0); }

if (argc==4)
   {
   if (!strcmp(argv[3],"-a")) option=1;
   else    { printf("\n ERROR:  option?  %s\n",argv[3]); exit(0);}
   }

/* compute the float transition probability matrix  */

for (k=j=0; j<DIM; j++) for (i=0; i<DIM; i++)
                          pam0[i][j]= (double) mdm[k++]/10000.;

for (i=0; i<DIM; i++) for (j=0; j<DIM; j++) pamn[i][j]=pam0[i][j];

for (n=N-1; n; n--)
    {
    Msqmul(pamn,pam0,store,DIM);
    for (i=0; i<DIM; i++) for (j=0; j<DIM; j++) pamn[i][j]=store[i][j];
    }

/* compute the pamnij sums  -> 1 ?     */

for (i=0; i<DIM; i++)
    { for (S=0., j=0; j<DIM; j++) S += pamn[i][j];
      printf (" pam %d  line  %d   sum: %lg\n",N, i, S);
    }

if (argc ==2)                    /* print out the nth pam */
   { for (i=0; i<DIM; i++)
       { for (j=0; j<DIM; j++) printf(" %4g", pamn[i][j]); putchar('\n');}
     exit(0);
   }



/* ************************************************************************* */
                   /* continue with transition matrix	*/

size=gt_Dmat(argv[2],T,&i,&j);
if (i != DIM || j !=DIM)
   { printf(" ERROR:  Transition matrix must be %d  x %d\n", DIM,DIM);
     exit(0); }
if (size < 0) { printf ("\n ERROR: not found %s\n", argv[2]); exit(0); }
if (size ==0) { printf ("\n ERROR: matrix format  %s\n",argv[2]); exit(0); }

/* compute sum transition matrix */
for (S=0.,i=0; i<DIM; i++)
            for (j=0; j<DIM; j++) S += T[i][j];
    printf ("\n Tij frequency sum: %lg\n", S);

if (!option)
   {
/* compute   T.M  -> store	*/
   Msqmul(T,pamn,store,DIM);

/* compute the T.Mij sum  -> 1 ?     */
   for (S=0.,i=0; i<DIM; i++)
       for (j=0; j<DIM; j++) S += store[i][j];
       printf ("\n TMij frequency sum: %lg\n", S);

/* compute     P(i,i)M(i,k)   ; P: frequency of i, M: pamn   */
   for (i=0; i<DIM; i++)
        for (k=0; k<DIM; k++) PM[i][k] = fq[i]*pamn[i][k];

/* compute PM sum    */
   for (S=0.,i=0; i<DIM; i++)
     for (j=0; j<DIM; j++) S += PM[i][j];
   printf ("\n PMij sum: %lg\n", S);

/* compute the final target frequency matrix  -> TM */
   Msqmul(PM,store,T,DIM);
  }
  else       /* ASYMETRICAL    [PAM]x [T]   */
  {

/* compute   MxT  -> store	*/
   Msqmul(pamn, T,store,DIM);

/* compute  final  target frequency matrix T:  P(i,i)MT(i,k) ; P: freq of i   */
   for (i=0; i<DIM; i++)
        for (k=0; k<DIM; k++) T[i][k] = fq[i]*store[i][k];
  }

/* compute the qij sum  -> 1 ?     */
   for (S=0.,i=0; i<DIM; i++)
     for (j=0; j<DIM; j++) S += T[i][j];
   printf ("\n Target frequency sum: %lg\n", S);

/* compute the pij sum  -> 1 ?     */
   for (S=0.,i=0; i<DIM; i++)
     for (j=0; j<DIM; j++) { PP[i][j]=fq[i]*fq[j]; S +=PP[i][j]; }
     printf ("\n pij frequency sum: %lg\n", S);

/* compute the relative entropy  */

   for (S=0.,i=0; i<DIM; i++)
     for (j=0; j<DIM; j++)
         S += T[i][j] * log(T[i][j]/PP[i][j]);
    printf("\n H = %lf\n",1.4427*S);                  /* 1/ln(2)    */

/* compute the log-odd score matrix */
for (i=0; i<DIM; i++)
    for (j=0; j<DIM; j++)
         T[i][j]=log(T[i][j]/(PP[i][j]));

/* print out the result  */

for (i=0; i<DIM; i++)
    { for (j=0; j<DIM; j++) printf(" %4.2f", T[i][j]); putchar('\n');}

for (lambda = 2.; lambda >0.001; lambda -= 0.1)
    { for (S=0.,i=0; i<DIM; i++)
          for (j=0; j<DIM; j++)
               S += PP[i][j]*exp(lambda*T[i][j]);
      printf ("\n lambda = %lf  (1-sum) = %lf\n", lambda, (1.-S));
    }

/* divide the scores by ln(2)/2  -> express as half-bit and -> DIM+2   */
for (i=0; i<DIM; i++)
    for (j=0; j<DIM; j++)
         MAT[i][j] = 2.88539*T[i][j];

/* compute additional lines and columns  */
/*  B = (D + N)/2  i=2,3         Z = (Q + E)/2   i=5,6   X = 0 */

for (j=0; j<DIM; j++)
    { MAT[DIM][j] = (T[2][j] + T[3][j])/2.	;	/* B  row */
      MAT[DIM+1][j] = (T[5][j] + T[6][j])/2.	;	/* Z  row */
      MAT[DIM+3][j] = 0.			;		/* X  row */
      MAT[j][DIM] = (T[j][2] + T[j][3])/2.	;	/* B  col */
      MAT[j][DIM+1] = (T[j][5] + T[j][5])/2.	;	/* B  col */
      MAT[j][DIM+2] = 0.			;		/* X  row */
    }
/* compute additional 3 x 3 matrix	*/
      MAT[DIM][DIM] = (T[2][2] + T[3][3])/2	;	   /* B diagonal */
      MAT[DIM+1][DIM+1] = (T[5][5] + T[6][6])/2	;	   /* Z diagonal */
      MAT[DIM+2][DIM+2] = 0.			;	   /* X diagonal */
      MAT[DIM+2][DIM] = 0.			;	   /* X, B */
      MAT[DIM+2][DIM+1] = 0.			;	   /* X  Z */
      MAT[DIM][DIM+2] = 0.			;	   /* B X */
      MAT[DIM+1][DIM+2] = 0.			;	   /* Z X */
      MAT[DIM][DIM+1] = (MAT[DIM][5] + MAT[DIM][6])/2.	;   /* B, Z      */
      MAT[DIM+1][DIM] = (MAT[DIM+1][2] + MAT[DIM+1][3])/2. ; /* Z, B      */

/* get the whole thing out     */

F=fopen("PAMxT.out", "w");
if (option)
  fprintf (F,"#    [PAM %d] [ %s ] scoring matrix\n", N, argv[2]);
else
  fprintf (F,"#    [PAM %d] [ %s ] [PAM %d]  scoring matrix\n", N, argv[2], N);
fprintf (F,"#\n");
fprintf 
(F,"  A  R  N  D  C  Q  E  G  H  I  L  K  M  F  P  S  T  W  Y  V  B  Z  X\n");
for (i=0; i<DIM+3; i++)
    { for (j=0; j<DIM+3; j++) fprintf(F," %2.0lf",MAT[i][j]);
      fputc((int)'\n',F);
    }
fclose(F);         

}

/* --------------------------------------------------------------------	*/
int Msqmul(mat1,mat2, mat3, I)
double mat1[][DIM], mat2[][DIM], mat3[][DIM];
int I;
{
int i,j,k;
for (i=0; i<I; i++)
  for (j=0; j<I; j++)
      for (k=0, mat3[i][j]=0.; k<I; k++) mat3[i][j] += mat1[i][k]*mat2[k][j];
}

/* -------------------------------------------------------------------	*/
/* int gt_Dmat()         Jean-Michel  Claverie   feb 1991   SGI		*/
/* Module to input a matrix  (float or int) into  a float array		*/
/* ARGUMENTS:								*/
/* IN:  fname      the matrix file name  , a closed string		*/
/* OUT:	M          double matrix, with max dimension set by MLIMIT	*/
/*      pI, pJ       *int dimensions : I lines,  J columns		*/
/* RETURNS:								*/
/*  -1 : for opening error; 0 for incorrect matrix; MLIMIT for overflow	*/
/* -------------------------------------------------------------------- */
/*             Matrix     format					*/
/* Leading comment lines are ignored, until the first DATA line.	*/
/* The first DATA line is defined as the first one with a valid scanf	*/
/* number as first word.						*/
/* The number of column is fixed from this 1st DATA line.		*/
/* Subsequent lines are presumed DATA lines, with their first J words   */
/* as matrix elements. Anything after the Jth word is a comment.	*/
/* Commentary trailing lines are not allowed.				*/
/* -------------------------------------------------------------------- */
/* ex:		comment line 1						*/
/* 		comment line 2						*/
/* 		comment line n						*/
/* 		n01  n02  n03  n04  comment				*/
/* 		n11  n12  n13  n14  comment				*/
/* 		nij ..............  comment				*/
/* -------------------------------------------------------------------- */	
int gt_Dmat(fname, M, pI, pJ)
char *fname;
int *pI,*pJ;
double M[][MLIMIT];
{
FILE *F;
char line[MAXLINE], *pline, word[100];
int i,j,next,sline;
float v;
if (NULL==(F=fopen(fname,"r")))
   { printf("\n ERROR: opening  %s\n",fname); return (-1); }

while (1){ if(-1 == fgt_line(F, MAXLINE,line)){ fclose(F); return 0; }
           if(sline == MAXLINE){ printf ("\n ERROR: line too long\n");
                                 fclose(F); return 0; }
           if(1==sscanf(line,"%lf", &v)) break;
         }
/*........ for the first line, only the J (unknown) first numbers are taken */
for (i=0,j=0, pline=line; (next=gt_1stword(pline,word)); pline += next, j++)
    { if(1!=sscanf(word,"%lf", &M[i][j])) break;   /* bail out, not a number */
      if (j==MLIMIT)           { fclose (F); return MLIMIT; }
    }
*pJ=j;
/*........ for the subsequent lines, only the first J words are taken */
for(i=1; (sline=fgt_line(F, MAXLINE, line)) !=-1 && sline < MAXLINE; i++)
   {for(j=0,pline=line;(next=gt_1stword(pline,word)) && j< *pJ;pline +=next,j++)
       { if(1!=sscanf(word,"%lf", &M[i][j])) break;
         if (j==MLIMIT)        { fclose (F); return MLIMIT; }
       }
    if (j != (*pJ))
       {printf("\n ERROR: column number ? %s\n",fname);fclose(F);return 0; }
    if (i==MLIMIT)     { fclose (F); return MLIMIT; }
   }
if(sline == MAXLINE){ printf ("\n ERROR: line too long\n"); fclose(F); return 0; }

*pI=i;    
return ((*pI > *pJ) ? *pI : *pJ );           /* return max dimension */
}


/* ---------------------------------------------------------------------*/
/* fgt_line(Fin, max,line)		JMC jan 1991       SGI		*/
/* #include <stdio.h> 							*/
/* function to input a line from an opened steam file into a char array,*/
/* with checking on the max length < max				*/
/* file "records" must be delimited by '\n'				*/
/* string is filled out until \n is encountered, or max-1 chars read	*/
/* or EOF occurs							*/
/* arguments:								*/
/*	Fin: FILE pointer to input stream				*/
/* 	max: array dimension [], max length including the final \0	*/
/*     line: string or char *						*/
/* return: the number of usable chars in the string (\0 not included)	*/
/*	 : or max if capacity was exceeded, line contains   		*/
/*	          max-1 usable characters.				*/
/*	 : 0 if NEW LINE was encountered first, line starts with '\0'	*/
/*	 :-1 if EOF was encountered, line starts with '\0'		*/
/* -------------------------------------------------------------------- */

fgt_line (Fin,max,line)
int max;
char line[];
FILE *Fin;

{
int c, i;

for (i=0; (c=fgetc(Fin))!= EOF && c!='\n' && i<(max-1);i++) line[i]=c;

line[i]='\0';          	/* replace NEW LINE or last by '\0', don't count */

if (c==EOF) return(-1); 			/* EOF encountered */
	
if (i==(max-1) && c!='\n') { while (fgetc(Fin)!='\n');     /* clean buffer */
			   return(max); }         	/* flag overflow   */
	return (i);
}  

/* sgt_line do the same things from a string */

/* -------------------------------------------------------------------	*/
/* MDtransp(c)     Jean-Michel  Claverie   feb 1991    SGI		*/
/* Module to transpose a double matrix within itself   			*/
/* ARGUMENTS:								*/
/* IN:    M, pI, pJ,  : M[I][J]   float matrix, dimension I, J		*/
/* OUT:	  transposed  M, NEW VALUES for calling I and J			*/
/* -------------------------------------------------------------------	*/
int MDtransp(M,pI,pJ)
double M[][MLIMIT];
int *pI,*pJ;
{
int i,j,It,Jt;
double v;
Jt= (*pI); It= (*pJ);
if (It>=Jt) 
   {	for (i=0 ; i< Jt; i++)
    	    for (j=i+1; j< It; j++)
                { v=M[j][i];
                  M[j][i]=M[i][j];
                  M[i][j]=v;
                }
   }
else
   {	for (j=0 ; j< It; j++)
    	    for (i=j+1; i< Jt; i++)
                { v=M[j][i];
                  M[j][i]=M[i][j];
                  M[i][j]=v;
                }
   }

*pI=It; *pJ=Jt;
}


/* -------------------------------------------------------------------	*/
/* int gt_1stword()         Jean-Michel Claverie      feb 1991    SGI	*/
/* module to read the first word of a closed string			*/
/* ARGUMENTS 								*/
/* IN:   line, a closed string						*/
/* OUT:  word, a closed string,  WITH sufficient space allocation	*/
/* RETURN: next: the relative position of the first characters NOT   	*/
/*         belonging to the first word in the string			*/
/*         thus, no more word in the string returns   0 (false)		*/
/* Note : the next call can then use : line + next, as a pointer to	*/
/*        initiate the search for the next 1st word in the string	*/
/* --------------------------------------------------------------------	*/
		
int gt_1stword (line,word)
char *line, *word;
{
char *p;
for ( p=line; isspace(*p) ; p++);  /* skeep leading blanks */
if (*p == '\0') return (0) ;		/* no word in line */

while (!isspace(*p) && *p!='\0'){ *word =(*p) ; word++; p++;}
*word='\0';
return (p-line);
}