LocalAlignmentStringKernel.cpp

Go to the documentation of this file.
00001 /*
00002  * Compute the local alignment kernel
00003  *
00004  * Largely based on LAkernel.c (version 0.3)
00005  *
00006  * Copyright 2003 Jean-Philippe Vert
00007  * Copyright 2005 Jean-Philippe Vert, Hiroto Saigo
00008  *
00009  * Shogun specific adjustments Written (W) 2007-2008 Soeren Sonnenburg
00010  * 
00011  * Reference:
00012  * H. Saigo, J.-P. Vert, T. Akutsu and N. Ueda, "Protein homology
00013  * detection using string alignment kernels", Bioinformatics,
00014  * vol.20, p.1682-1689, 2004.
00015  * 
00016  * This program is free software; you can redistribute it and/or modify
00017  * it under the terms of the GNU General Public License as published by
00018  * the Free Software Foundation; either version 3 of the License, or
00019  * (at your option) any later version.
00020  */
00021 
00022 #include <stdlib.h>
00023 #include <stdio.h>
00024 #include <math.h>
00025 #include <ctype.h>
00026 #include <string.h>
00027 #include "kernel/LocalAlignmentStringKernel.h"
00028 
00029 using namespace shogun;
00030 
00031 /****************/
00032 /* The alphabet */
00033 /****************/
00034 
00035 #define NAA 20                                  /* Number of amino-acids */
00036 #define NLET 26                                 /* Number of letters in the alphabet */
00037 const char* CLocalAlignmentStringKernel::aaList= "ARNDCQEGHILKMFPSTWYV";    /* The list of amino acids */
00038 
00039 /*****************/
00040 /* SW parameters */
00041 /*****************/
00042 
00043 #define OPENING 12                              /* Gap opening penalty */
00044 #define EXTENSION 2                             /* Gap extension penalty */
00045 
00046 /* mutation matrix */
00047 const int32_t CLocalAlignmentStringKernel::blosum[] = {
00048   6,
00049  -2,   8,
00050  -2,  -1,   9,
00051   -3,  -2,   2,   9,
00052   -1,  -5,  -4,  -5,  13,
00053   -1,   1,   0,   0,  -4,   8,
00054   -1,   0,   0,   2,  -5,   3,   7,
00055   0,  -3,  -1,  -2,  -4,  -3,  -3,   8,
00056   -2,   0,   1,  -2,  -4,   1,   0,  -3,  11,
00057   -2,  -5,  -5,  -5,  -2,  -4,  -5,  -6,  -5,   6,
00058   -2,  -3,  -5,  -5,  -2,  -3,  -4,  -5,  -4,   2,   6,
00059   -1,   3,   0,  -1,  -5,   2,   1,  -2,  -1,  -4,  -4,   7,
00060   -1,  -2,  -3,  -5,  -2,  -1,  -3,  -4,  -2,   2,   3,  -2,   8,
00061   -3,  -4,  -5,  -5,  -4,  -5,  -5,  -5,  -2,   0,   1,  -5,   0,   9,
00062   -1,  -3,  -3,  -2,  -4,  -2,  -2,  -3,  -3,  -4,  -4,  -2,  -4,  -5,  11,
00063   2,  -1,   1,   0,  -1,   0,   0,   0,  -1,  -4,  -4,   0,  -2,  -4,  -1,   6,
00064   0,  -2,   0,  -2,  -1,  -1,  -1,  -2,  -3,  -1,  -2,  -1,  -1,  -3,  -2,   2,   7,
00065   -4,  -4,  -6,  -6,  -3,  -3,  -4,  -4,  -4,  -4,  -2,  -4,  -2,   1,  -6,  -4,  -4,  16,
00066   -3,  -3,  -3,  -5,  -4,  -2,  -3,  -5,   3,  -2,  -2,  -3,  -1,   4,  -4,  -3,  -2,   3,  10,
00067   0,  -4,  -4,  -5,  -1,  -3,  -4,  -5,  -5,   4,   1,  -3,   1,  -1,  -4,  -2,   0,  -4,  -2,   6};
00068 
00069 /* Index corresponding to the (i,j) entry (i,j=0..19) in the blosum matrix */
00070 #define BINDEX(i,j) (((i)>(j))?(j)+(((i)*(i+1))/2):(i)+(((j)*(j+1))/2))
00071 
00072 /*********************
00073  * Kernel parameters *
00074  *********************/
00075 
00076 #define SCALING 0.1           /* Factor to scale all SW parameters */
00077 
00078 /* If you want to compute the sum over all local alignments (to get a valid kernel), uncomment the following line : */
00079 /* If x=log(a) and y=log(b), compute log(a+b) : */
00080 /*
00081 #define LOGP(x,y) (((x)>(y))?(x)+log1p(exp((y)-(x))):(y)+log1p(exp((x)-(y))))
00082 */
00083 
00084 #define LOGP(x,y) LogSum(x,y)
00085 
00086 /* OR if you want to compute the score of the best local alignment (to get the SW score by Viterbi), uncomment the following line : */
00087 /*
00088 #define LOGP(x,y) (((x)>(y))?(x):(y))
00089 */
00090 
00091 /* Usefule constants */
00092 #define LOG0 -10000          /* log(0) */
00093 #define INTSCALE 1000.0      /* critical for speed and precise computation*/
00094 
00095 int32_t CLocalAlignmentStringKernel::logsum_lookup[LOGSUM_TBL];
00096 
00097 CLocalAlignmentStringKernel::CLocalAlignmentStringKernel(int32_t size)
00098 : CStringKernel<char>(size), initialized(false)
00099 {
00100     scaled_blosum=new int32_t[sizeof(blosum)];
00101     init_logsum();
00102     initialize();
00103 }
00104 
00105 CLocalAlignmentStringKernel::CLocalAlignmentStringKernel(
00106     CStringFeatures<char>* l, CStringFeatures<char>* r)
00107 : CStringKernel<char>(10), initialized(false)
00108 {
00109     scaled_blosum=new int32_t[sizeof(blosum)];
00110     init_logsum();
00111     initialize();
00112     init(l, r);
00113 }
00114 
00115 CLocalAlignmentStringKernel::~CLocalAlignmentStringKernel()
00116 {
00117     cleanup();
00118 }
00119 
00120 bool CLocalAlignmentStringKernel::init(CFeatures* l, CFeatures* r)
00121 {
00122     CStringKernel<char>::init(l, r);
00123     initialized = true;
00124     return init_normalizer();
00125 }
00126 
00127 void CLocalAlignmentStringKernel::cleanup()
00128 {
00129     delete[] scaled_blosum;
00130     scaled_blosum=NULL;
00131 
00132     free(isAA);
00133     isAA=NULL;
00134     free(aaIndex);
00135     aaIndex=NULL;
00136 
00137     CKernel::cleanup();
00138 }
00139 
00140 /* LogSum - default log funciotion. fast, but not exact */
00141 /* LogSum2 - precise, but slow. Note that these two functions need different figure types  */
00142 
00143 void CLocalAlignmentStringKernel::init_logsum(void){
00144   int32_t i;
00145   for (i = 0; i < LOGSUM_TBL; i++) 
00146     logsum_lookup[i] = (int32_t) (INTSCALE*
00147                    (log(1.+exp( (float32_t) -i/INTSCALE))));
00148 }
00149 
00150 int32_t CLocalAlignmentStringKernel::LogSum(int32_t p1, int32_t p2)
00151 {
00152     int32_t diff;
00153     static int32_t firsttime=1;
00154 
00155     if (firsttime)
00156     {
00157         init_logsum();
00158         firsttime =0;
00159     }
00160 
00161     diff=p1-p2;
00162     if (diff>=LOGSUM_TBL) return p1;
00163     else if (diff<=-LOGSUM_TBL) return p2;
00164     else if (diff>0) return p1+logsum_lookup[diff];
00165     else return p2+logsum_lookup[-diff];
00166 }
00167 
00168 
00169 float32_t CLocalAlignmentStringKernel::LogSum2(float32_t p1, float32_t p2)
00170 {
00171   if (p1 > p2)
00172     return (p1-p2 > 50.) ? p1 : p1 + log(1. + exp(p2-p1));
00173   else
00174     return (p2-p1 > 50.) ? p2 : p2 + log(1. + exp(p1-p2));
00175 }
00176 
00177 
00178 void CLocalAlignmentStringKernel::initialize(void)
00179      /* Initialize all static variables. This function should be called once before computing the first pair HMM score */
00180 {
00181   register int32_t i;
00182 
00183   /* Initialization of the array which gives the position of each amino-acid in the set of amino-acid */
00184   if ((aaIndex=(int32_t *)calloc(NLET,sizeof(int32_t))) == NULL)
00185     SG_ERROR("run out o memory");
00186   for (i=0;i<NAA;i++) 
00187     aaIndex[aaList[i]-'A']=i;
00188   
00189   /* Initialization of the array which indicates whether a char is an amino-acid */
00190   if ((isAA=(int32_t *)calloc(256,sizeof(int32_t))) == NULL)
00191     SG_ERROR("run out of memory");
00192   for (i=0;i<NAA;i++) 
00193     isAA[(int32_t)aaList[i]]=1;
00194 
00195   /* Scale the blossum matrix */
00196   for (i=0 ; i<NAA*(NAA+1)/2; i++)
00197       scaled_blosum[i] = (int32_t) floor(blosum[i]*SCALING*INTSCALE);
00198 
00199 
00200   /* Scale of gap penalties */
00201   opening = (int32_t) floor(OPENING * SCALING*INTSCALE);
00202   extension = (int32_t) floor(EXTENSION * SCALING*INTSCALE);
00203 }
00204 
00205 
00206 
00207 /* Implementation of the
00208  * convolution kernel which generalizes the Smith-Waterman algorithm
00209  */
00210 float64_t CLocalAlignmentStringKernel::LAkernelcompute(
00211     int32_t* aaX, int32_t* aaY, /* the two amino-acid sequences (as sequences of indexes in [0..NAA-1] indicating the position of the amino-acid in the variable 'aaList') */
00212     int32_t nX, int32_t nY /* the lengths of both sequences */)
00213 {
00214    register int32_t
00215     i,j,                /* loop indexes */
00216     cur, old,           /* to indicate the array to use (0 or 1) */
00217     curpos, frompos;    /* position in an array */
00218 
00219    int32_t
00220     *logX,           /* arrays to store the log-values of each state */
00221     *logY,
00222     *logM,
00223     *logX2,
00224     *logY2,
00225 
00226     aux , aux2;/* , aux3 , aux4 , aux5;*/
00227   int32_t
00228     cl;                /* length of a column for the dynamic programming */
00229 
00230   /*
00231   printf("now computing pairHMM between %d and %d:\n",nX,nY);
00232   for (i=0;i<nX;printf("%d ",aaX[i++]));
00233   printf("\n and \n");
00234   for (i=0;i<nY;printf("%d ",aaY[i++]));
00235   printf("\n");
00236   */
00237 
00238   /* Initialization of the arrays */
00239   /* Each array stores two successive columns of the (nX+1)x(nY+1) table used in dynamic programming */
00240   cl = nY+1;           /* each column stores the positions in the aaY sequence, plus a position at zero */
00241 
00242   logM=new int32_t[2*cl];
00243   logX=new int32_t[2*cl];
00244   logY=new int32_t[2*cl];
00245   logX2=new int32_t[2*cl];
00246   logY2=new int32_t[2*cl];
00247 
00248   /************************************************/
00249   /* First iteration : initialization of column 0 */
00250   /************************************************/
00251   /* The log=proabilities of each state are initialized for the first column (x=0,y=0..nY) */
00252 
00253   for (j=0;j<cl;j++) {
00254     logM[j]=LOG0;
00255     logX[j]=LOG0;
00256     logY[j]=LOG0;
00257     logX2[j]=LOG0;
00258     logY2[j]=LOG0;
00259 
00260   }
00261 
00262   /* Update column order */
00263   cur = 1;      /* Indexes [0..cl-1] are used to process the next column */
00264   old = 0;      /* Indexes [cl..2*cl-1] were used for column 0 */
00265 
00266 
00267   /************************************************/
00268   /* Next iterations : processing columns 1 .. nX */
00269   /************************************************/
00270 
00271   /* Main loop to vary the position in aaX : i=1..nX */
00272   for (i=1;i<=nX;i++) {
00273 
00274     /* Special update for positions (i=1..nX,j=0) */
00275     curpos = cur*cl;                  /* index of the state (i,0) */
00276     logM[curpos] = LOG0; 
00277     logX[curpos] = LOG0; 
00278     logY[curpos] = LOG0; 
00279     logX2[curpos] = LOG0; 
00280     logY2[curpos] = LOG0; 
00281 
00282     /* Secondary loop to vary the position in aaY : j=1..nY */
00283     for (j=1;j<=nY;j++) {
00284 
00285       curpos = cur*cl + j;            /* index of the state (i,j) */
00286 
00287       /* Update for states which emit X only */
00288       /***************************************/
00289 
00290       frompos = old*cl + j;            /* index of the state (i-1,j) */
00291       
00292       /* State RX */
00293       logX[curpos] = LOGP( - opening + logM[frompos] , - extension + logX[frompos] );
00294       /*      printf("%.5f\n",logX[curpos]);*/
00295       /*      printf("%.5f\n",logX_B[curpos]);*/
00296       /* State RX2 */
00297       logX2[curpos] = LOGP( logM[frompos] , logX2[frompos] );
00298 
00299       /* Update for states which emit Y only */
00300       /***************************************/
00301 
00302       frompos = cur*cl + j-1;          /* index of the state (i,j-1) */
00303 
00304       /* State RY */
00305       aux = LOGP( - opening + logM[frompos] , - extension + logY[frompos] );
00306       logY[curpos] = LOGP( aux , - opening + logX[frompos] );
00307 
00308       /* State RY2 */
00309       aux = LOGP( logM[frompos] , logY2[frompos] );
00310       logY2[curpos] = LOGP( aux , logX2[frompos] );
00311 
00312       /* Update for states which emit X and Y */
00313       /****************************************/
00314 
00315       frompos = old*cl + j-1;          /* index of the state (i-1,j-1) */
00316 
00317       aux = LOGP( logX[frompos] , logY[frompos] );
00318       aux2 = LOGP( 0 , logM[frompos] );
00319       logM[curpos] = LOGP( aux , aux2 ) + scaled_blosum[ BINDEX( aaX[i-1] , aaY[j-1] ) ];
00320       
00321       /*      printf("i=%d , j=%d\nM=%.5f\nX=%.5f\nY=%.5f\nX2=%.5f\nY2=%.5f\n",i,j,logM[curpos],logX[curpos],logY[curpos],logX2[curpos],logY2[curpos]);
00322        */
00323 
00324     }  /* end of j=1:nY loop */
00325 
00326 
00327     /* Update the culumn order */
00328     cur = 1-cur;
00329     old = 1-old;
00330 
00331   }  /* end of j=1:nX loop */
00332 
00333 
00334   /* Termination */
00335   /***************/
00336 
00337   curpos = old*cl + nY;                /* index of the state (nX,nY) */
00338   aux = LOGP( logX2[curpos] , logY2[curpos] );
00339   aux2 = LOGP( 0 , logM[curpos] );
00340   /*  kernel_value = LOGP( aux , aux2 );*/
00341 
00342   /* Memory release */
00343     delete[] logM;
00344     delete[] logX;
00345     delete[] logY;
00346     delete[] logX2;
00347     delete[] logY2;
00348 
00349   /* Return the logarithm of the kernel */
00350   return (float32_t)LOGP(aux,aux2)/INTSCALE;
00351 }
00352 
00353 /********************/
00354 /* Public functions */
00355 /********************/
00356 
00357 
00358 /* Return the log-probability of two sequences x and y under a pair HMM model */
00359 /* x and y are strings of aminoacid letters, e.g., "AABRS" */
00360 float64_t CLocalAlignmentStringKernel::compute(int32_t idx_x, int32_t idx_y)
00361 {
00362   int32_t *aax,*aay;  /* to convert x and y into sequences of amino-acid indexes */
00363   int32_t lx=0,ly=0;       /* lengths of x and y */
00364   int32_t i,j;
00365 
00366   /* If necessary, initialize static variables */
00367   if (isAA == NULL)
00368     initialize();
00369 
00370   bool free_x, free_y;
00371   char* x=((CStringFeatures<char>*) lhs)->get_feature_vector(idx_x, lx, free_x);
00372   char* y=((CStringFeatures<char>*) rhs)->get_feature_vector(idx_y, ly, free_y);
00373   ASSERT(x && y);
00374 
00375   if ((lx<1) || (ly<1))
00376     SG_ERROR("empty chain");
00377 
00378   /* Create aax and aay */
00379 
00380   if ((aax=(int32_t *)calloc(lx,sizeof(int32_t))) == NULL)
00381     SG_ERROR("run out of memory");
00382   if ((aay=(int32_t *)calloc(ly,sizeof(int32_t))) == NULL)
00383     SG_ERROR("run out of memory");
00384 
00385   /* Extract the characters corresponding to aminoacids and keep their indexes */
00386 
00387   j=0;
00388   for (i=0 ; i<lx ; i++) 
00389     if (isAA[toupper(x[i])])
00390       aax[j++] = aaIndex[toupper(x[i])-'A'];
00391   lx = j;
00392   j=0;
00393   for (i=0 ; i<ly ; i++)
00394     if (isAA[toupper(y[i])])
00395       aay[j++] = aaIndex[toupper(y[i])-'A'];
00396   ly = j;
00397 
00398 
00399   /* Compute the pair HMM score */
00400   float64_t result=LAkernelcompute(aax,aay,lx,ly);
00401 
00402   /* Release memory */
00403   free(aax);
00404   free(aay);
00405 
00406   ((CStringFeatures<char>*) lhs)->free_feature_vector(x, idx_x, free_x);
00407   ((CStringFeatures<char>*) rhs)->free_feature_vector(y, idx_y, free_y);
00408 
00409   return result;
00410 }

SHOGUN Machine Learning Toolbox - Documentation