binaries/src/fasta34/tatstats.c

   1 #include <stdlib.h>
   2 #include <stdio.h>
   3 #include <math.h>
   4 #include <string.h>
   5
   6 #include "defs.h"
   7 #include "param.h"
   8 #include "tatstats.h"
   9
  10 #ifndef PCOMPLIB
  11 #include "mw.h"
  12 #else
  13 #include "p_mw.h"
  14 #endif
  15
  16 /* calc_priors() - calculate frequencies of amino-acids, possibly with counts */
  17 /* generate_tatprobs() - build the table of score probabilities if the
  18    sequences are not too long */
  19
  20 double
  21 det(double a11, double a12, double a13,
  22     double a21, double a22, double a23,
  23     double a31, double a32, double a33);
  24
  25 double power(double r, int p)
  26 {
  27   double tr;
  28   int neg;
  29
  30   if (r==0.0) return(p==0?1.0:0.0);
  31   if (neg = p<0) p = -p;
  32   tr = 1.0;
  33   while (p>0) {
  34     if (p & 1) tr *= r;
  35     p >>= 1;
  36     if (p) r *= r;
  37   }
  38   return((neg? 1.0/tr: tr));
  39 }
  40
  41 double
  42 factorial (int a, int b) {
  43
  44   double res = 1.0;
  45
  46   if(a == 0) { return 1.0; }
  47
  48   while(a > b) {
  49     res *= (double) a;
  50     a--;
  51   }
  52
  53   return res;
  54 }
  55
  56 void
  57 calc_priors(double *priors,
  58             struct pstruct *ppst,
  59             struct f_struct *f_str,
  60             const unsigned char *aa1, int n1,
  61             int pseudocts)
  62 {
  63   long counts[25], sum;
  64   int i;
  65
  66   if(n1 == 0 && f_str->priors[1] > 0.0) {
  67     for(i = 1 ; i <= ppst->nsq ; i++) {
  68       priors[i] = f_str->priors[i];
  69     }
  70     return;
  71   }
  72
  73   if(n1 == 0) {
  74     if (ppst->dnaseq==SEQT_PROT ) {
  75
  76       /* Robinson & Robinson residue counts from Stephen Altschul */
  77       counts[ 1] = 35155;    /*   A  */
  78       counts[ 2] = 23105;    /*   R  */
  79       counts[ 3] = 20212;    /*   N  */
  80       counts[ 4] = 24161;    /*   D  */
  81       counts[ 5] =  8669;    /*   C  */
  82       counts[ 6] = 19208;    /*   Q  */
  83       counts[ 7] = 28354;    /*   E  */
  84       counts[ 8] = 33229;    /*   G  */
  85       counts[ 9] =  9906;    /*   H  */
  86       counts[10] = 23161;    /*   I  */
  87       counts[11] = 40625;    /*   L  */
  88       counts[12] = 25872;    /*   K  */
  89       counts[13] = 10101;    /*   M  */
  90       counts[14] = 17367;    /*   F  */
  91       counts[15] = 23435;    /*   P  */
  92       counts[16] = 32070;    /*   S  */
  93       counts[17] = 26311;    /*   T  */
  94       counts[18] =  5990;    /*   W  */
  95       counts[19] = 14488;    /*   Y  */
  96       counts[20] = 29012;    /*   V  */
  97       counts[21] =     0;    /*   B  */
  98       counts[22] =     0;    /*   Z  */
  99       counts[23] =     0;    /*   X   */
 100       counts[24] =     0;    /*   *   */
 101     }
 102     else { /* SEQT_DNA */
 103       counts[1] = 250;
 104       counts[2] = 250;
 105       counts[3] = 250;
 106       counts[4] = 250;
 107       for (i=5; i<=ppst->nsq; i++) counts[i]=0;
 108     }
 109   } else {
 110     memset(&counts[0], 0, sizeof(counts));
 111
 112     for(i = 0 ; i < n1 ; i++) {
 113       if(aa1[i] > ppst->nsq || aa1[i] < 1) continue;
 114       counts[aa1[i]]++;
 115     }
 116   }
 117
 118   sum = 0;
 119   for(i = 1 ; i <= ppst->nsq ; i++) sum += counts[i];
 120
 121   for(i = 1 ; i <= ppst->nsq ; i++) {
 122     if(n1 == 0) {
 123       priors[i] = (double) counts[i] / (double) sum;
 124     } else {
 125       priors[i] = ( ((double) pseudocts * f_str->priors[i]) + (double) counts[i] ) / ( (double) sum + (double) pseudocts );
 126     }
 127   }
 128
 129   return;
 130 }
 131
 132 int
 133 max_score(int *scores, int nsq) {
 134
 135   int max, i;
 136
 137   max = -BIGNUM;
 138   for ( i = 1 ; i <= nsq ; i++ ) {
 139     if (scores[i] > max) max = scores[i];
 140   }
 141
 142   return max;
 143 }
 144
 145 int
 146 min_score(int *scores, int nsq) {
 147
 148   int min, i;
 149
 150   min = BIGNUM;
 151   for (i = 1 ; i <= nsq ; i++ ) {
 152     if (scores[i] < min) min = scores[i];
 153   }
 154
 155   return min;
 156 }
 157
 158 double
 159 calc_tatusov ( struct slink *last,
 160                struct slink *this,
 161                const unsigned char *aa0, int n0,
 162                const unsigned char *aa1, int n1,
 163                int **pam2, int nsq,
 164                struct f_struct *f_str,
 165                int pseudocts,
 166                int do_opt,
 167                int zsflag
 168                )
 169 {
 170   int i, is, j, k;
 171
 172   double *priors, my_priors[MAXSQ], tatprob, left_tatprob, right_tatprob;
 173   unsigned char *query = NULL;
 174   int length, maxlength, sumlength, sumscore, tmp, seg;
 175   int start, stop;
 176   struct slink *sl;
 177   int N;
 178   double *tatprobsptr;
 179
 180 #if defined(FASTS) || defined(FASTM)
 181   int index = 0;
 182   int notokay = 0;
 183 #endif
 184
 185   struct tat_str *oldtat = NULL, *newtat = NULL;
 186
 187 #if defined(FASTS) || defined(FASTM)
 188   start = this->vp->start - this->vp->dp + f_str->noff;
 189   stop = this->vp->stop - this->vp->dp + f_str->noff;
 190   tmp = stop - start + 1;
 191 #else
 192   /*
 193     FASTF alignments can also hang off the end of library sequences,
 194     but no query residues are used up in the process, but we have to
 195     keep track of which are
 196   */
 197   tmp = 0;
 198   for(i = 0, j = 0 ; i < n0 ; i++) {
 199     if (this->vp->used[i] == 1) {tmp++; }
 200   }
 201 #endif
 202
 203   sumlength = maxlength = length = tmp;
 204   seg = 1;
 205   sumscore = this->vp->score;
 206
 207 #if defined(FASTS) || defined(FASTM)
 208   if(f_str->aa0b[start] == start && f_str->aa0e[stop] == stop) {
 209     index |= (1 << f_str->aa0i[start]);
 210   } else {
 211     notokay |= (1 << f_str->aa0i[start]);
 212   }
 213 #endif
 214
 215   for(sl = last; sl != NULL ; sl = sl->prev) {
 216
 217 #if defined(FASTS) || defined(FASTM)
 218     start = sl->vp->start - sl->vp->dp + f_str->noff;
 219     stop = sl->vp->stop - sl->vp->dp + f_str->noff;
 220     tmp = stop - start + 1;
 221 #else
 222     tmp = 0;
 223     for(i = 0, j = 0 ; i < n0 ; i++) {
 224       if(sl->vp->used[i] == 1) {
 225         tmp++;
 226       }
 227     }
 228 #endif
 229     sumlength += tmp;
 230     maxlength = tmp > maxlength ? tmp : maxlength;
 231     seg++;
 232     sumscore += sl->vp->score;
 233
 234 #if defined(FASTS) || defined(FASTM)
 235     if(f_str->aa0b[start] == start && f_str->aa0e[stop] == stop) {
 236       index |= (1 << f_str->aa0i[start]);
 237     } else {
 238       notokay |= (1 << f_str->aa0i[start]);
 239     }
 240 #endif
 241
 242   }
 243
 244   tatprob = -1.0;
 245
 246 #if defined(FASTS) || defined(FASTM)
 247
 248   /* for T?FASTS, we try to use what we've precalculated: */
 249
 250   /* with z = 3, do_opt is true, but we can use precalculated - with
 251      all other z's we can use precalculated only if !do_opt */
 252   if(!notokay && f_str->tatprobs != NULL) {
 253     /* create our own newtat and copy f_str's tat into it */
 254     index--;
 255
 256     newtat = (struct tat_str *) malloc(sizeof(struct tat_str));
 257     if(newtat == NULL) {
 258       fprintf(stderr, "Couldn't calloc memory for newtat.\n");
 259       exit(1);
 260     }
 261
 262     memcpy(newtat, f_str->tatprobs[index], sizeof(struct tat_str));
 263
 264     newtat->probs = (double *) calloc(f_str->tatprobs[index]->highscore - f_str->tatprobs[index]->lowscore + 1, sizeof(double));
 265     if(newtat->probs == NULL) {
 266       fprintf(stderr, "Coudln't calloc memory for newtat->probs.\n");
 267       exit(1);
 268     }
 269
 270     memcpy(newtat->probs, f_str->tatprobs[index]->probs,
 271            (f_str->tatprobs[index]->highscore - f_str->tatprobs[index]->lowscore + 1) * sizeof(double));
 272
 273
 274     tatprob = f_str->intprobs[index][sumscore - f_str->tatprobs[index]->lowscore];
 275
 276   } else { /* we need to recalculate from scratch */
 277 #endif
 278
 279     /* for T?FASTF, we're always recalculating from scratch: */
 280
 281     query = (unsigned char *) calloc(length, sizeof(unsigned char));
 282     if(query == NULL) {
 283       fprintf(stderr, "Couldn't calloc memory for query.\n");
 284       exit(1);
 285     }
 286
 287 #if defined(FASTS) || defined(FASTM)
 288     start = this->vp->start - this->vp->dp + f_str->noff;
 289     for(i = 0, j = 0 ; i < length ; i++) {
 290       query[j++] = aa0[start + i];
 291     }
 292 #else
 293     for(i = 0, j = 0 ; i < n0 ; i++) {
 294       if (this->vp->used[i] == 1) {query[j++] = aa0[i];}
 295     }
 296 #endif
 297
 298     /*  calc_priors - not currently implemented for aa1 dependent */
 299     /*
 300     if( (do_opt && zsflag == 2) || zsflag == 4 ) {
 301       priors = &my_priors[0];
 302       calc_priors(priors, f_str, aa1, n1, pseudocts);
 303     } else {
 304       priors = f_str->priors;
 305     }
 306     */
 307
 308     priors = f_str->priors;
 309     oldtat = (last != NULL ? last->tat : NULL);
 310
 311     generate_tatprobs(query, 0, length - 1, priors, pam2, nsq, &newtat, oldtat);
 312
 313     free(query);
 314 #if defined(FASTS) || defined(FASTM)
 315   } /* close the FASTS-specific if-else from above */
 316 #endif
 317
 318   this->newtat = newtat;
 319
 320   if(tatprob < 0.0) { /* hasn't been set by precalculated FASTS intprobs */
 321
 322     /* integrate probabilities >= sumscore */
 323     tatprobsptr = newtat->probs;
 324
 325     is = i = newtat->highscore - newtat->lowscore;
 326     N = sumscore - newtat->lowscore;
 327
 328     right_tatprob = 0;
 329     for ( ;  i >= N; i--) {
 330       right_tatprob += tatprobsptr[i];
 331     }
 332
 333     left_tatprob = tatprobsptr[0];
 334     for (i = 1 ; i < N ; i++ ) {
 335       left_tatprob += tatprobsptr[i];
 336     }
 337
 338     if (right_tatprob < left_tatprob) {tatprob = right_tatprob;}
 339     else {tatprob = 1.0 - left_tatprob;}
 340
 341     tatprob /= (right_tatprob+left_tatprob);
 342   }
 343
 344   if (maxlength > 0) {
 345     n1 += 2 * (maxlength - 1);
 346   }
 347
 348 #ifndef FASTM
 349   tatprob *= factorial(n1 - sumlength + seg, n1 - sumlength);
 350 #else
 351   tatprob *= power(n1 - sumlength,seg)/(1<<seg);
 352 #endif
 353
 354   if(tatprob > 0.01)
 355     tatprob = 1.0 - exp(-tatprob);
 356
 357   return tatprob;
 358 }
 359
 360 void
 361 generate_tatprobs(const unsigned char *query,
 362                   int begin,
 363                   int end,
 364                   double *priors,
 365                   int **pam2,
 366                   int nsq,
 367                   struct tat_str **tatarg,
 368                   struct tat_str *oldtat)
 369 {
 370
 371   int i, j, k, l, m, n, N, highscore, lowscore;
 372   int *lowrange = NULL, *highrange = NULL;
 373   double *probs = NULL, *newprobs = NULL, *priorptr, tmp;
 374   struct tat_str *tatprobs = NULL;
 375   int *pamptr, *pamptrsave;
 376
 377   if((tatprobs = (struct tat_str *) calloc(1, sizeof(struct tat_str)))==NULL) {
 378     fprintf(stderr, "Couldn't allocate individual tatprob struct.\n");
 379     exit(1);
 380   }
 381
 382   n = end - begin + 1;
 383
 384   if ( (lowrange = (int *) calloc(n, sizeof(int))) == NULL ) {
 385     fprintf(stderr, "Couldn't allocate memory for lowrange.\n");
 386     exit(1);
 387   }
 388
 389   if ( (highrange = (int *) calloc(n, sizeof(int))) == NULL ) {
 390     fprintf(stderr, "Couldn't allocate memory for highrange.\n");
 391     exit(1);
 392   }
 393
 394   /* calculate the absolute highest and lowest score possible for this */
 395   /* segment.  Also, set the range we need to iterate over at each position */
 396   /* in the query: */
 397   if(oldtat == NULL) {
 398     highscore = lowscore = 0;
 399   } else {
 400     highscore = oldtat->highscore;
 401     lowscore = oldtat->lowscore;
 402   }
 403
 404   for ( i = 0 ; i < n ; i++ ) {
 405
 406     if (query[begin+i] == 0) break;
 407
 408     highscore =
 409       (highrange[i] = highscore + max_score(pam2[query[begin + i]], nsq));
 410
 411     lowscore =
 412       (lowrange[i] = lowscore + min_score(pam2[query[begin + i]], nsq));
 413
 414     /*
 415     fprintf(stderr, "i: %d, max: %d, min: %d, high[i]: %d, low[i]: %d, high: %d, low: %d, char: %d\n",
 416             i,
 417             max_score(pam2[query[begin + i]], nsq),
 418             min_score(pam2[query[begin + i]], nsq),
 419             highrange[i], lowrange[i],
 420             highscore, lowscore, query[begin + i]);
 421     */
 422   }
 423
 424   /* allocate an array of probabilities for all possible scores */
 425   /* i.e. if highest score possible is 50 and lowest score possible */
 426   /* is -20, then there are 50 - (-20) + 1 = 71 possible different */
 427   /* scores (including 0): */
 428   N = highscore - lowscore;
 429   if ( (probs = (double *) calloc(N + 1, sizeof(double))) == NULL ) {
 430     fprintf(stderr, "Couldn't allocate probability matrix : %d.\n", N + 1);
 431     exit(1);
 432   }
 433
 434   if(oldtat == NULL) {
 435     /* for the first position, iterate over the only possible scores, */
 436     /* summing the priors for the amino acids that can yield each score. */
 437     pamptr = pam2[query[begin]];
 438     for ( i = 1 ; i <= nsq ; i++ ) {
 439       if(priors[i] > 0.0) {
 440         probs[(pamptr[i] - lowscore)] += priors[i];
 441       }
 442     }
 443   } else {
 444     /* Need to copy the data out of oldtat->probs into probs */
 445     memcpy( &probs[oldtat->lowscore - lowscore],
 446             oldtat->probs,
 447             (oldtat->highscore - oldtat->lowscore + 1) * sizeof(double));
 448   }
 449
 450   if ( (newprobs = (double *) calloc(N + 1, sizeof(double))) == NULL ) {
 451     fprintf(stderr, "Couldn't allocate newprobs matrix.\n");
 452     exit(1);
 453   }
 454
 455   /* now for each remaining residue in the segment ... */
 456   for ( i = (oldtat == NULL ? 1 : 0) ; i < n ; i++ ) {
 457
 458     pamptrsave = pam2[query[begin + i]];
 459
 460     /* ... calculate new probability distribution .... */
 461
 462     /* ... for each possible score (limited to current range) ... */
 463     for ( j = lowrange[i] - lowscore,
 464             k = highrange[i] - lowscore ;
 465           j <= k ;
 466           j++ ) {
 467
 468       tmp = 0.0;
 469       pamptr = &pamptrsave[1];
 470       priorptr = &priors[1];
 471       /* ... for each of the possible alignment scores at this position ... */
 472       for ( l = 1 ;
 473             l <= nsq ;
 474             l++) {
 475
 476         /* make sure we don't go past highest possible score, or past
 477            the lowest possible score; not sure why this can happen */
 478         m = j - *pamptr++;
 479         if ( m <= N && m >= 0 ) {
 480           /* update the probability of getting score j: */
 481           tmp += probs[m] * *priorptr++;
 482         }
 483       }
 484       newprobs[j] += tmp;
 485     }
 486
 487     /* save the new set of probabilities, get rid of old; we don't
 488        necessarily have to copy/clear all N+1 slots, we could use
 489        high/low score boundaries -- not sure that's worth the
 490        effort. */
 491     memcpy(probs, newprobs, (N + 1) * sizeof(double));
 492     memset(newprobs, 0, (N + 1) * sizeof(double));
 493   }
 494
 495   free(newprobs);
 496   free(highrange);
 497   free(lowrange);
 498
 499   tatprobs->probs = probs;
 500   /*  tatprobs->intprobs = intprobs; */
 501   tatprobs->lowscore = lowscore;
 502   tatprobs->highscore = highscore;
 503
 504   *tatarg = tatprobs;
 505 }
 506