binaries/src/fasta34/print_pssm.c

   1 /* print_pssm.c - 21-Jan-2005
   2
   3    copyright (c) 2005 - William R. Pearson and the University of Virginia
   4
   5    read a binary PSSM checkpoint file from blastpgp, and produce an ascii
   6    formatted file
   7 */
   8
   9 #include <stdio.h>
  10 #include <stdlib.h>
  11 #include <sys/types.h>
  12 #include <math.h>
  13 #include <string.h>
  14
  15 #include "defs.h"
  16 #include "mm_file.h"
  17 #include "param.h"
  18
  19 #include "uascii.h"
  20 #include "upam.h"
  21
  22 void initenv(int, char **, struct pstruct *, char *);
  23 void read_pssm();
  24 void alloc_pam();
  25 int **alloc_pam2p();
  26 void initpam2();
  27 void fill_pam();
  28 double get_lambda();
  29
  30 extern int optind;
  31 extern char *optarg;
  32
  33 main(int argc, char **argv) {
  34
  35   char *aa0;
  36   char libstr[MAX_FN];
  37   char qname[MAX_FN];
  38   int sq0off;
  39   int i, n0;
  40   FILE *fp;
  41   struct pstruct pst, *ppst;
  42
  43   /* stuff from initfa.c/h_init() */
  44
  45   memcpy(qascii,aascii,sizeof(qascii));
  46
  47   /* initialize a pam matrix */
  48   ppst = &pst;
  49   strncpy(ppst->pamfile,"BL50",MAX_FN);
  50   standard_pam(ppst->pamfile,ppst,0,0);
  51
  52   /* this is always protein by default */
  53   ppst->nsq = naa;
  54   ppst->nsqx = naax;
  55   for (i=0; i<=ppst->nsqx; i++) {
  56     ppst->sq[i] = aa[i];
  57     ppst->hsq[i] = haa[i];
  58     ppst->sqx[i]=aax[i];        /* sq = aa */
  59     ppst->hsqx[i]=haax[i];      /* hsq = haa */
  60   }
  61   ppst->sq[ppst->nsqx+1] = ppst->sqx[ppst->nsqx+1] = '\0';
  62
  63   if ((aa0 = calloc(MAXTST,sizeof(char)))==NULL) {
  64     fprintf(stderr,"Cannot allocate aa0\n");
  65     exit(1);
  66   }
  67
  68   initenv(argc, argv, &pst, qname);
  69   alloc_pam(pst.nsq+1,pst.nsq+1, &pst);
  70   initpam2(&pst);
  71
  72   n0 = getseq (qname, qascii, aa0, MAXTST, libstr,&sq0off);
  73
  74   if (!pst.pam_pssm) {
  75     fprintf(stderr," ** ERROR ** No -P PSSM provided\n");
  76   }
  77   else {
  78     ppst->pam2p[0] = alloc_pam2p(n0,pst.nsq);
  79     ppst->pam2p[1] = alloc_pam2p(n0,pst.nsq);
  80     if ((fp = fopen(pst.pgpfile,"rb"))!=NULL) {
  81       read_pssm(aa0, n0, pst.nsq, pst.pamscale,fp,ppst);
  82     }
  83   }
  84 }
  85
  86 void
  87 initenv(int argc, char **argv, struct pstruct *ppst, char *qname) {
  88   char copt;
  89
  90   pascii = aascii;
  91
  92   while ((copt = getopt(argc, argv, "P:s:"))!=EOF) {
  93     switch (copt) {
  94       case 'P':
  95         strncpy(ppst->pgpfile,optarg,MAX_FN);
  96         ppst->pgpfile[MAX_FN-1]='\0';
  97         ppst->pam_pssm = 1;
  98         break;
  99
 100       case 's':
 101         strncpy (ppst->pamfile, optarg, 120);
 102         ppst->pamfile[120-1]='\0';
 103         if (!standard_pam(ppst->pamfile,ppst,0, 0)) {
 104           initpam (ppst->pamfile, ppst);
 105         }
 106         ppst->pam_set=1;
 107         break;
 108     }
 109   }
 110   optind--;
 111
 112   if (argc - optind > 1) strncpy(qname, argv[optind+1], MAX_FN);
 113 }
 114
 115
 116 /*
 117    *aa0 - query sequence
 118    n0   - length
 119    pamscale - scaling for pam matrix - provided by apam.c, either
 120               0.346574 = ln(2)/2 (P120, BL62) or
 121               0.231049 = ln(2)/3 (P250, BL50)
 122 */
 123
 124 #define N_EFFECT 20
 125
 126 void
 127 read_pssm(unsigned char *aa0, int n0, int nsq, double pamscale, FILE *fp, struct pstruct *ppst) {
 128   int i, j, len;
 129   int qi, rj;
 130   int **pam2p;
 131   int first, too_high;
 132   char *query;
 133   double freq, **freq2d, lambda, new_lambda;
 134   double scale, scale_high, scale_low;
 135
 136   pam2p = ppst->pam2p[0];
 137
 138   if(1 != fread(&len, sizeof(int), 1, fp)) {
 139     fprintf(stderr, "error reading from checkpoint file: %d\n", len);
 140     exit(1);
 141   }
 142
 143   if(len != n0) {
 144     fprintf(stderr, "profile length (%d) and query length (%d) don't match!\n",
 145             len,n0);
 146     exit(1);
 147   }
 148
 149   /* read over query sequence stored in BLAST profile */
 150   if(NULL == (query = (char *) calloc(len, sizeof(char)))) {
 151     fprintf(stderr, "Couldn't allocate memory for query!\n");
 152     exit(1);
 153   }
 154
 155   if(len != fread(query, sizeof(char), len, fp)) {
 156     fprintf(stderr, "Couldn't read query sequence from profile: %s\n", query);
 157     exit(1);
 158   }
 159
 160   printf("%d\n%s\n",len,query);
 161
 162   /* currently we don't do anything with query; ideally, we should
 163      check to see that it actually matches aa0 ... */
 164
 165   /* quick 2d array alloc: */
 166   if((freq2d = (double **) calloc(n0, sizeof(double *))) == NULL) {
 167     fprintf(stderr, "Couldn't allocate memory for frequencies!\n");
 168     exit(1);
 169   }
 170
 171   if((freq2d[0] = (double *) calloc(n0 * N_EFFECT, sizeof(double))) == NULL) {
 172     fprintf(stderr, "Couldn't allocate memory for frequencies!\n");
 173     exit(1);
 174   }
 175
 176   /* a little pointer arithmetic to fill out 2d array: */
 177   for (qi = 1 ; qi < n0 ; qi++) {
 178     freq2d[qi] = freq2d[0] + (N_EFFECT * qi);
 179   }
 180
 181   for (qi = 0 ; qi < n0 ; qi++) {
 182     printf("%c",query[qi]);
 183     for (rj = 0 ; rj < N_EFFECT ; rj++) {
 184       if(1 != fread(&freq, sizeof(double), 1, fp)) {
 185         fprintf(stderr, "Error while reading frequencies!\n");
 186         exit(1);
 187       }
 188       printf(" %8.7g",freq*10.0);
 189
 190       if (freq > 1e-12) {
 191         freq = log(freq /((double) (rrcounts[rj+1])/(double) rrtotal));
 192         freq /= pamscale; /* this gets us close to originial pam scores */
 193         freq2d[qi][rj] = freq;
 194       }
 195       else {freq2d[qi][rj] = freq;}
 196     }
 197     printf("\n");
 198   }
 199
 200
 201   /* now figure out the right scale */
 202   scale = 1.0;
 203   lambda = get_lambda(ppst->pam2[0], 20, 20, "\0ARNDCQEGHILKMFPSTWYV");
 204
 205   /* should be near 1.0 because of our initial scaling by ppst->pamscale */
 206   fprintf(stderr, "real_lambda: %g\n", lambda);
 207
 208   /* get initial high/low scale values: */
 209   first = 1;
 210   while (1) {
 211     fill_pam(pam2p, n0, 20, freq2d, scale);
 212     new_lambda = get_lambda(pam2p, n0, 20, query);
 213
 214     if (new_lambda > lambda) {
 215       if (first) {
 216         first = 0;
 217         scale = scale_high = 1.0 + 0.05;
 218         scale_low = 1.0;
 219         too_high = 1;
 220       } else {
 221         if (!too_high) break;
 222         scale = (scale_high += scale_high - 1.0);
 223       }
 224     } else if (new_lambda > 0) {
 225       if (first) {
 226         first = 0;
 227         scale_high = 1.0;
 228         scale = scale_low = 1.0 - 0.05;
 229         too_high = 0;
 230       } else {
 231         if (too_high) break;
 232         scale = (scale_low += scale_low - 1.0);
 233       }
 234     } else {
 235       fprintf(stderr, "new_lambda (%g) <= 0; matrix has positive average score", new_lambda);
 236       exit(1);
 237     }
 238   }
 239
 240   /* now do binary search between low and high */
 241   for (i = 0 ; i < 10 ; i++) {
 242     scale = 0.5 * (scale_high + scale_low);
 243     fill_pam(pam2p, n0, 20, freq2d, scale);
 244     new_lambda = get_lambda(pam2p, n0, 20, query);
 245
 246     if (new_lambda > lambda) scale_low = scale;
 247     else scale_high = scale;
 248   }
 249
 250   scale = 0.5 * (scale_high + scale_low);
 251   fill_pam(pam2p, n0, 20, freq2d, scale);
 252
 253   fprintf(stderr, "final scale: %g\n", scale);
 254
 255   for (qi = 0 ; qi < n0 ; qi++) {
 256     fprintf(stderr, "%4d %c:  ", qi+1, query[qi]);
 257     for (rj = 1 ; rj <= 20 ; rj++) {
 258       fprintf(stderr, "%4d", pam2p[qi][rj]);
 259     }
 260     fprintf(stderr, "\n");
 261   }
 262
 263   free(freq2d[0]);
 264   free(freq2d);
 265
 266   free(query);
 267 }
 268
 269 /*
 270  * alloc_pam(): allocates memory for the 2D pam matrix as well
 271  * as for the integer array used to transmit the pam matrix
 272  */
 273 void
 274 alloc_pam (int d1, int d2, struct pstruct *ppst)
 275 {
 276   int     i, *d2p;
 277
 278   if ((ppst->pam2[0] = (int **) malloc (d1 * sizeof (int *))) == NULL) {
 279      fprintf(stderr,"Cannot allocate 2D pam matrix: %d",d1);
 280      exit(1);
 281   }
 282
 283   if ((ppst->pam2[1] = (int **) malloc (d1 * sizeof (int *))) == NULL) {
 284      fprintf(stderr,"Cannot allocate 2D pam matrix: %d",d1);
 285      exit(1);
 286   }
 287
 288   if ((d2p = pam12 = (int *) malloc (d1 * d2 * sizeof (int))) == NULL) {
 289      fprintf(stderr,"Cannot allocate 2D pam matrix: %d",d1);
 290      exit(1);
 291    }
 292
 293    for (i = 0; i < d1; i++, d2p += d2)
 294       ppst->pam2[0][i] = d2p;
 295
 296    if ((d2p=pam12x= (int *) malloc (d1 * d2 * sizeof (int))) == NULL) {
 297      fprintf(stderr,"Cannot allocate 2d pam matrix: %d",d2);
 298      exit(1);
 299    }
 300
 301    for (i = 0;  i < d1; i++, d2p += d2)
 302       ppst->pam2[1][i] = d2p;
 303 }
 304
 305 void
 306 fill_pam(int **pam2p, int n0, int nsq, double **freq2d, double scale) {
 307   int i, j;
 308   double freq;
 309
 310   /* fprintf(stderr, "scale: %g\n", scale); */
 311
 312   /* now fill in the pam matrix: */
 313   for (i = 0 ; i < n0 ; i++) {
 314     for (j = 1 ; j <=nsq ; j++) {
 315       freq = scale * freq2d[i][j-1];
 316       if ( freq < 0.0) freq -= 0.5;
 317       else freq += 0.5;
 318       pam2p[i][j] = (int)(freq);
 319     }
 320   }
 321 }
 322
 323 /*
 324  *  initpam2(struct pstruct pst): Converts 1-D pam matrix to 2-D
 325  */
 326 void initpam2 (struct pstruct *ppst)
 327 {
 328    int i, j, k, nsq, pam_xx, pam_xm;
 329    int sa_x, sa_t, tmp;
 330
 331    nsq = ppst->nsq;
 332    sa_x = pascii['X'];
 333    sa_t = pascii['*'];
 334
 335    ppst->pam2[0][0][0] = -BIGNUM;
 336    ppst->pam_h = -1; ppst->pam_l = 1;
 337
 338    k = 0;
 339    for (i = 1; i <= nsq; i++) {
 340      ppst->pam2[0][0][i] = ppst->pam2[0][i][0] = -BIGNUM;
 341      for (j = 1; j <= i; j++) {
 342        ppst->pam2[0][j][i] = ppst->pam2[0][i][j] = pam[k++] - ppst->pamoff;
 343        if (ppst->pam_l > ppst->pam2[0][i][j]) ppst->pam_l =ppst->pam2[0][i][j];
 344        if (ppst->pam_h < ppst->pam2[0][i][j]) ppst->pam_h =ppst->pam2[0][i][j];
 345      }
 346    }
 347
 348    ppst->nt_align = (ppst->dnaseq== SEQT_DNA || ppst->dnaseq == SEQT_RNA);
 349
 350    if (ppst->dnaseq == SEQT_RNA) {
 351      tmp = ppst->pam2[0][nascii['G']][nascii['G']] - 1;
 352      ppst->pam2[0][nascii['A']][nascii['G']] =
 353        ppst->pam2[0][nascii['C']][nascii['T']] =
 354        ppst->pam2[0][nascii['C']][nascii['U']] = tmp;
 355    }
 356
 357    if (ppst->pam_x_set) {
 358      for (i=1; i<=nsq; i++) {
 359        ppst->pam2[0][sa_x][i] = ppst->pam2[0][i][sa_x]=ppst->pam_xm;
 360        ppst->pam2[0][sa_t][i] = ppst->pam2[0][i][sa_t]=ppst->pam_xm;
 361      }
 362      ppst->pam2[0][sa_x][sa_x]=ppst->pam_xx;
 363      ppst->pam2[0][sa_t][sa_t]=ppst->pam_xm;
 364    }
 365    else {
 366      ppst->pam_xx = ppst->pam2[0][sa_x][sa_x];
 367      ppst->pam_xm = ppst->pam2[0][1][sa_x];
 368    }
 369 }
 370
 371 double
 372 get_lambda(int **pam2p, int n0, int nsq, char *aa0) {
 373   double lambda, H;
 374   double *pr, tot, sum;
 375   int i, ioff, j, min, max;
 376
 377   /* get min and max scores */
 378   min = BIGNUM;
 379   max = -BIGNUM;
 380   if(pam2p[0][1] == -BIGNUM) {
 381     ioff = 1;
 382     n0++;
 383   } else {
 384     ioff = 0;
 385   }
 386
 387   for (i = ioff ; i < n0 ; i++) {
 388     for (j = 1; j <= nsq ; j++) {
 389       if (min > pam2p[i][j])
 390         min = pam2p[i][j];
 391       if (max < pam2p[i][j])
 392         max = pam2p[i][j];
 393     }
 394   }
 395
 396   /*  fprintf(stderr, "min: %d\tmax:%d\n", min, max); */
 397
 398   if ((pr = (double *) calloc(max - min + 1, sizeof(double))) == NULL) {
 399     fprintf(stderr, "Couldn't allocate memory for score probabilities: %d\n", max - min + 1);
 400     exit(1);
 401   }
 402
 403   tot = (double) rrtotal * (double) rrtotal * (double) n0;
 404   for (i = ioff ; i < n0 ; i++) {
 405     for (j = 1; j <= nsq ; j++) {
 406       pr[pam2p[i][j] - min] +=
 407         (double) ((double) rrcounts[aascii[aa0[i]]] * (double) rrcounts[j]) / tot;
 408     }
 409   }
 410
 411   sum = 0.0;
 412   for(i = 0 ; i <= max-min ; i++) {
 413     sum += pr[i];
 414     /*    fprintf(stderr, "%3d: %g %g\n", i+min, pr[i], sum); */
 415   }
 416   /*  fprintf(stderr, "sum: %g\n", sum); */
 417
 418   for(i = 0 ; i <= max-min ; i++) { pr[i] /= sum; }
 419
 420   if (!karlin(min, max, pr, &lambda, &H)) {
 421     fprintf(stderr, "Karlin lambda estimation failed\n");
 422   }
 423
 424   /*   fprintf(stderr, "lambda: %g\n", lambda); */
 425   free(pr);
 426
 427   return lambda;
 428 }
 429
 430 int **
 431 alloc_pam2p(int len, int nsq) {
 432   int i;
 433   int **pam2p;
 434
 435   if ((pam2p = (int **)calloc(len,sizeof(int *)))==NULL) {
 436     fprintf(stderr," Cannot allocate pam2p: %d\n",len);
 437     return NULL;
 438   }
 439
 440   if((pam2p[0] = (int *)calloc((nsq+1)*len,sizeof(int)))==NULL) {
 441     fprintf(stderr, "Cannot allocate pam2p[0]: %d\n", (nsq+1)*len);
 442     free(pam2p);
 443     return NULL;
 444   }
 445
 446   for (i=1; i<len; i++) {
 447     pam2p[i] = pam2p[0] + (i*(nsq+1));
 448   }
 449
 450   return pam2p;
 451 }
 452
 453 void free_pam2p(int **pam2p) {
 454   if (pam2p) {
 455     free(pam2p[0]);
 456     free(pam2p);
 457   }
 458 }
 459