/*-------------------------------------------------------------

Copyright (C) 2000 Peter Clote, Sebastian Will. 
All Rights Reserved.

written by Peter Clote
modified by Sebastian Will

Permission to use, copy, modify, and distribute this
software and its documentation for NON-COMMERCIAL purposes
and without fee is hereby granted provided that this
copyright notice appears in all copies.


THE AUTHOR AND PUBLISHER MAKE NO REPRESENTATIONS OR
WARRANTIES ABOUT THE SUITABILITY OF THE SOFTWARE, EITHER
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE, OR NON-INFRINGEMENT. THE AUTHORS
AND PUBLISHER SHALL NOT BE LIABLE FOR ANY DAMAGES SUFFERED
BY LICENSEE AS A RESULT OF USING, MODIFYING OR DISTRIBUTING
THIS SOFTWARE OR ITS DERIVATIVES.

-------------------------------------------------------------*/

/*---------------------------------------------------
jensenShannon.c

P. Clote 19 Dec 1999

This program opens a file, computes length N of file and
prints out the Jensen Shannon divergence when the
file is segmented between 1,..., i*step and
i*step+1,...,N.

High divergence indicates that the left and right
segments are more homogeneous with respect to themselves
than with respect to the whole.

This could be combined with an entropy plot, dinucleotide
entropy plot, or pair potential plot.

----------------------------------------------------*/


#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#include<time.h>

#define STEP 100
   // print JensenShannon divergence for multiplicities of STEP
#define LOG2(x)  (log(x)/log(2))



/*------- prototypes --------*/
char nucleotide(char ch);
int isPurine(char ch);
int isPyrimidine(char ch);
double computeEntropy(int pur, int pyr, int size);


int main(int argc, char *argv[]){
   FILE *in;
   time_t time_seed = time(0);
   int purine,pyrimidine; 
   int purineU,pyrimidineU,purineV,pyrimidineV; 
   int i,j,num,n,nU,nV;
      // entire genome W of length n is segmented into U,V
   int ch;
   double HW,HU,HV,divergence; // entropy values


   /*********  ERROR CHECKING in command line arguments *********/
   if (argc != 2) {
       fprintf(stderr,"%s filename\n",argv[0]);
       exit(1);
   }
   if ((in = fopen(argv[1],"r")) == NULL) {
       fprintf(stderr,"%s not a readable file.\n", argv[1]);
       exit(1);
   }

   srand(time_seed); 
   /*--------------------------------------------
   Read file twice through, and each time make same
   random assignments to A,C,G,T for incomplete data.
   --------------------------------------------*/


   n=nU=nV=purine=pyrimidine=purineU=pyrimidineU=purineV=pyrimidineV=0;
     // initialize counts to 0

   while ((ch = fgetc(in)) != EOF) {
     char nucl=nucleotide(ch);
     if ( isPurine(nucl) )
         purine++;
     else if (isPyrimidine(nucl))
         pyrimidine++;
     if (nucl) n++;
   }   
   rewind(in);
   srand(time_seed); 
      /*--------------------------------------------
      Reset to same random seed to make same
      random assignments to A,C,G,T for incomplete data.
      --------------------------------------------*/


   HW = computeEntropy(purine,pyrimidine,n);
   printf("purine: %d pyrimidine: %d genome length: %d\n",
           purine,pyrimidine,n);
   printf("Entropy H(W): %lf\n",HW);


   // Compute JensenShannon divergence at segments
   // U = 1,...,i*STEP and V = i*STEP+1,....,n
   purineV=purine;
   pyrimidineV=pyrimidine;
   nV=n;

   num = (n-1)/STEP;
   for (i=0;i<num;i++) {
     for(j=0; j<STEP; j++){
       char nucl;
       do nucl = nucleotide(ch = fgetc(in));
       while (nucl==0 && ch!=EOF);
       if (ch==EOF) printf("e\n");

       if ( isPurine(nucl) ) {
           purineU++;
           purineV--;
           nU++;
           nV--;
           }
       else { // nucl must be pyrimidine
           pyrimidineU++;
           pyrimidineV--;
           nU++;
           nV--;
           }
     }// end for j

     HU = computeEntropy(purineU,pyrimidineU,nU);
     HV = computeEntropy(purineV,pyrimidineV,nV);
     divergence = HW - (double) nU/n * HU - (double) nV/n * HV;
     // printf("HU: %lf\t",HU);
     // printf("HV: %lf\n",HV);
     printf("%lf\n",divergence);
     
     }//end for 

    fclose(in);
} //end main


int isPurine(char ch){
   return( ch == 'A' || ch == 'G' );
}

int isPyrimidine(char ch){
   return( ch == 'C' || ch == 'T' );
}


/* nucleotide(ch)
   function to handle incomplete data in sequences
   if base ch is a completely determined base A,C,G, or T this base 
   is returned.
   If ch is a symbol for a uncertain base, a base is returned
   at random according to the symbol
   If ch is not a valid symbol at all, 0 is returned
 */

char nucleotide(char ch){
   int x;

   switch (toupper(ch)) {
      case 'A' : return(ch); 
      case 'C' : return(ch); 
      case 'G' : return(ch); 
      case 'T' : return(ch); 

      case 'N' : x = rand()%4;
                 if ( x == 0 ) 
                    return('A');   
                 else if ( x == 1)
                    return('C');  
                 else if ( x == 2)
                    return('G');  
                 else 
                    return('T');  

      case 'R' : if ( rand()%2 )  // purine
                    return('A');   
                 else
                    return('G');
      case 'Y' : if ( rand()%2 )  // pyrimidine
                    return('C');   
                 else
                    return('T');
      case 'W' : if ( rand()%2 )  // weak
                    return('A');   
                 else
                    return('T');
      case 'S' : if ( rand()%2 )  // strong
                    return('C');   
                 else
                    return('G');
      case 'M' : if ( rand()%2 ) 
                    return('A');   
                 else
                    return('C');  
      case 'K' : if ( rand()%2 ) 
                    return('G');   
                 else
                    return('T');  
      case 'V' : x = rand()%3;
                 if ( x == 0 ) 
                    return('A');   
                 else if ( x == 1)
                    return('C');  
                 else 
                    return('G');  
      case 'H' : x = rand()%3;
                 if ( x == 0 ) 
                    return('A');   
                 else if ( x == 1)
                    return('C');  
                 else 
                    return('T');  
      case 'D' : x = rand()%3;
                 if ( x == 0 ) 
                    return('A');   
                 else if ( x == 1)
                    return('G');  
                 else 
                    return('T');  
      case 'B' : x = rand()%3;
                 if ( x == 0 ) 
                    return('C');   
                 else if ( x == 1)
                    return('G');  
                 else 
                    return('T');
   default: return 0;
     } // end switch

}//end nucleotide



double computeEntropy(int pur, int pyr, int size){
   double p,q;
   if((pur<=0)||(pyr<=0)||(size<=0))
       return (0);

   p = (double) pur/size;
   q = (double) pyr/size;
   return (-p*LOG2(p) - q*LOG2(q));
}