1 /* print_pssm.c - 21-Jan-2005
3 copyright (c) 2005 - William R. Pearson and the University of Virginia
5 read a binary PSSM checkpoint file from blastpgp, and produce an ascii
11 #include <sys/types.h>
22 void initenv(int, char **, struct pstruct *, char *);
33 main(int argc, char **argv) {
41 struct pstruct pst, *ppst;
43 /* stuff from initfa.c/h_init() */
45 memcpy(qascii,aascii,sizeof(qascii));
47 /* initialize a pam matrix */
49 strncpy(ppst->pamfile,"BL50",MAX_FN);
50 standard_pam(ppst->pamfile,ppst,0,0);
52 /* this is always protein by default */
55 for (i=0; i<=ppst->nsqx; i++) {
57 ppst->hsq[i] = haa[i];
58 ppst->sqx[i]=aax[i]; /* sq = aa */
59 ppst->hsqx[i]=haax[i]; /* hsq = haa */
61 ppst->sq[ppst->nsqx+1] = ppst->sqx[ppst->nsqx+1] = '\0';
63 if ((aa0 = calloc(MAXTST,sizeof(char)))==NULL) {
64 fprintf(stderr,"Cannot allocate aa0\n");
68 initenv(argc, argv, &pst, qname);
69 alloc_pam(pst.nsq+1,pst.nsq+1, &pst);
72 n0 = getseq (qname, qascii, aa0, MAXTST, libstr,&sq0off);
75 fprintf(stderr," ** ERROR ** No -P PSSM provided\n");
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);
87 initenv(int argc, char **argv, struct pstruct *ppst, char *qname) {
92 while ((copt = getopt(argc, argv, "P:s:"))!=EOF) {
95 strncpy(ppst->pgpfile,optarg,MAX_FN);
96 ppst->pgpfile[MAX_FN-1]='\0';
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);
112 if (argc - optind > 1) strncpy(qname, argv[optind+1], MAX_FN);
117 *aa0 - query sequence
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)
127 read_pssm(unsigned char *aa0, int n0, int nsq, double pamscale, FILE *fp, struct pstruct *ppst) {
133 double freq, **freq2d, lambda, new_lambda;
134 double scale, scale_high, scale_low;
136 pam2p = ppst->pam2p[0];
138 if(1 != fread(&len, sizeof(int), 1, fp)) {
139 fprintf(stderr, "error reading from checkpoint file: %d\n", len);
144 fprintf(stderr, "profile length (%d) and query length (%d) don't match!\n",
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");
155 if(len != fread(query, sizeof(char), len, fp)) {
156 fprintf(stderr, "Couldn't read query sequence from profile: %s\n", query);
160 printf("%d\n%s\n",len,query);
162 /* currently we don't do anything with query; ideally, we should
163 check to see that it actually matches aa0 ... */
165 /* quick 2d array alloc: */
166 if((freq2d = (double **) calloc(n0, sizeof(double *))) == NULL) {
167 fprintf(stderr, "Couldn't allocate memory for frequencies!\n");
171 if((freq2d[0] = (double *) calloc(n0 * N_EFFECT, sizeof(double))) == NULL) {
172 fprintf(stderr, "Couldn't allocate memory for frequencies!\n");
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);
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");
188 printf(" %8.7g",freq*10.0);
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;
195 else {freq2d[qi][rj] = freq;}
201 /* now figure out the right scale */
203 lambda = get_lambda(ppst->pam2[0], 20, 20, "\0ARNDCQEGHILKMFPSTWYV");
205 /* should be near 1.0 because of our initial scaling by ppst->pamscale */
206 fprintf(stderr, "real_lambda: %g\n", lambda);
208 /* get initial high/low scale values: */
211 fill_pam(pam2p, n0, 20, freq2d, scale);
212 new_lambda = get_lambda(pam2p, n0, 20, query);
214 if (new_lambda > lambda) {
217 scale = scale_high = 1.0 + 0.05;
221 if (!too_high) break;
222 scale = (scale_high += scale_high - 1.0);
224 } else if (new_lambda > 0) {
228 scale = scale_low = 1.0 - 0.05;
232 scale = (scale_low += scale_low - 1.0);
235 fprintf(stderr, "new_lambda (%g) <= 0; matrix has positive average score", new_lambda);
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);
246 if (new_lambda > lambda) scale_low = scale;
247 else scale_high = scale;
250 scale = 0.5 * (scale_high + scale_low);
251 fill_pam(pam2p, n0, 20, freq2d, scale);
253 fprintf(stderr, "final scale: %g\n", scale);
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]);
260 fprintf(stderr, "\n");
270 * alloc_pam(): allocates memory for the 2D pam matrix as well
271 * as for the integer array used to transmit the pam matrix
274 alloc_pam (int d1, int d2, struct pstruct *ppst)
278 if ((ppst->pam2[0] = (int **) malloc (d1 * sizeof (int *))) == NULL) {
279 fprintf(stderr,"Cannot allocate 2D pam matrix: %d",d1);
283 if ((ppst->pam2[1] = (int **) malloc (d1 * sizeof (int *))) == NULL) {
284 fprintf(stderr,"Cannot allocate 2D pam matrix: %d",d1);
288 if ((d2p = pam12 = (int *) malloc (d1 * d2 * sizeof (int))) == NULL) {
289 fprintf(stderr,"Cannot allocate 2D pam matrix: %d",d1);
293 for (i = 0; i < d1; i++, d2p += d2)
294 ppst->pam2[0][i] = d2p;
296 if ((d2p=pam12x= (int *) malloc (d1 * d2 * sizeof (int))) == NULL) {
297 fprintf(stderr,"Cannot allocate 2d pam matrix: %d",d2);
301 for (i = 0; i < d1; i++, d2p += d2)
302 ppst->pam2[1][i] = d2p;
306 fill_pam(int **pam2p, int n0, int nsq, double **freq2d, double scale) {
310 /* fprintf(stderr, "scale: %g\n", scale); */
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;
318 pam2p[i][j] = (int)(freq);
324 * initpam2(struct pstruct pst): Converts 1-D pam matrix to 2-D
326 void initpam2 (struct pstruct *ppst)
328 int i, j, k, nsq, pam_xx, pam_xm;
335 ppst->pam2[0][0][0] = -BIGNUM;
336 ppst->pam_h = -1; ppst->pam_l = 1;
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];
348 ppst->nt_align = (ppst->dnaseq== SEQT_DNA || ppst->dnaseq == SEQT_RNA);
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;
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;
362 ppst->pam2[0][sa_x][sa_x]=ppst->pam_xx;
363 ppst->pam2[0][sa_t][sa_t]=ppst->pam_xm;
366 ppst->pam_xx = ppst->pam2[0][sa_x][sa_x];
367 ppst->pam_xm = ppst->pam2[0][1][sa_x];
372 get_lambda(int **pam2p, int n0, int nsq, char *aa0) {
374 double *pr, tot, sum;
375 int i, ioff, j, min, max;
377 /* get min and max scores */
380 if(pam2p[0][1] == -BIGNUM) {
387 for (i = ioff ; i < n0 ; i++) {
388 for (j = 1; j <= nsq ; j++) {
389 if (min > pam2p[i][j])
391 if (max < pam2p[i][j])
396 /* fprintf(stderr, "min: %d\tmax:%d\n", min, max); */
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);
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;
412 for(i = 0 ; i <= max-min ; i++) {
414 /* fprintf(stderr, "%3d: %g %g\n", i+min, pr[i], sum); */
416 /* fprintf(stderr, "sum: %g\n", sum); */
418 for(i = 0 ; i <= max-min ; i++) { pr[i] /= sum; }
420 if (!karlin(min, max, pr, &lambda, &H)) {
421 fprintf(stderr, "Karlin lambda estimation failed\n");
424 /* fprintf(stderr, "lambda: %g\n", lambda); */
431 alloc_pam2p(int len, int nsq) {
435 if ((pam2p = (int **)calloc(len,sizeof(int *)))==NULL) {
436 fprintf(stderr," Cannot allocate pam2p: %d\n",len);
440 if((pam2p[0] = (int *)calloc((nsq+1)*len,sizeof(int)))==NULL) {
441 fprintf(stderr, "Cannot allocate pam2p[0]: %d\n", (nsq+1)*len);
446 for (i=1; i<len; i++) {
447 pam2p[i] = pam2p[0] + (i*(nsq+1));
453 void free_pam2p(int **pam2p) {