--- /dev/null
+/* print_pssm.c - 21-Jan-2005
+
+ copyright (c) 2005 - William R. Pearson and the University of Virginia
+
+ read a binary PSSM checkpoint file from blastpgp, and produce an ascii
+ formatted file
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <sys/types.h>
+#include <math.h>
+#include <string.h>
+
+#include "defs.h"
+#include "mm_file.h"
+#include "param.h"
+
+#include "uascii.h"
+#include "upam.h"
+
+void initenv(int, char **, struct pstruct *, char *);
+void read_pssm();
+void alloc_pam();
+int **alloc_pam2p();
+void initpam2();
+void fill_pam();
+double get_lambda();
+
+extern int optind;
+extern char *optarg;
+
+main(int argc, char **argv) {
+
+ char *aa0;
+ char libstr[MAX_FN];
+ char qname[MAX_FN];
+ int sq0off;
+ int i, n0;
+ FILE *fp;
+ struct pstruct pst, *ppst;
+
+ /* stuff from initfa.c/h_init() */
+
+ memcpy(qascii,aascii,sizeof(qascii));
+
+ /* initialize a pam matrix */
+ ppst = &pst;
+ strncpy(ppst->pamfile,"BL50",MAX_FN);
+ standard_pam(ppst->pamfile,ppst,0,0);
+
+ /* this is always protein by default */
+ ppst->nsq = naa;
+ ppst->nsqx = naax;
+ for (i=0; i<=ppst->nsqx; i++) {
+ ppst->sq[i] = aa[i];
+ ppst->hsq[i] = haa[i];
+ ppst->sqx[i]=aax[i]; /* sq = aa */
+ ppst->hsqx[i]=haax[i]; /* hsq = haa */
+ }
+ ppst->sq[ppst->nsqx+1] = ppst->sqx[ppst->nsqx+1] = '\0';
+
+ if ((aa0 = calloc(MAXTST,sizeof(char)))==NULL) {
+ fprintf(stderr,"Cannot allocate aa0\n");
+ exit(1);
+ }
+
+ initenv(argc, argv, &pst, qname);
+ alloc_pam(pst.nsq+1,pst.nsq+1, &pst);
+ initpam2(&pst);
+
+ n0 = getseq (qname, qascii, aa0, MAXTST, libstr,&sq0off);
+
+ if (!pst.pam_pssm) {
+ fprintf(stderr," ** ERROR ** No -P PSSM provided\n");
+ }
+ else {
+ ppst->pam2p[0] = alloc_pam2p(n0,pst.nsq);
+ ppst->pam2p[1] = alloc_pam2p(n0,pst.nsq);
+ if ((fp = fopen(pst.pgpfile,"rb"))!=NULL) {
+ read_pssm(aa0, n0, pst.nsq, pst.pamscale,fp,ppst);
+ }
+ }
+}
+
+void
+initenv(int argc, char **argv, struct pstruct *ppst, char *qname) {
+ char copt;
+
+ pascii = aascii;
+
+ while ((copt = getopt(argc, argv, "P:s:"))!=EOF) {
+ switch (copt) {
+ case 'P':
+ strncpy(ppst->pgpfile,optarg,MAX_FN);
+ ppst->pgpfile[MAX_FN-1]='\0';
+ ppst->pam_pssm = 1;
+ break;
+
+ case 's':
+ strncpy (ppst->pamfile, optarg, 120);
+ ppst->pamfile[120-1]='\0';
+ if (!standard_pam(ppst->pamfile,ppst,0, 0)) {
+ initpam (ppst->pamfile, ppst);
+ }
+ ppst->pam_set=1;
+ break;
+ }
+ }
+ optind--;
+
+ if (argc - optind > 1) strncpy(qname, argv[optind+1], MAX_FN);
+}
+
+
+/*
+ *aa0 - query sequence
+ n0 - length
+ pamscale - scaling for pam matrix - provided by apam.c, either
+ 0.346574 = ln(2)/2 (P120, BL62) or
+ 0.231049 = ln(2)/3 (P250, BL50)
+*/
+
+#define N_EFFECT 20
+
+void
+read_pssm(unsigned char *aa0, int n0, int nsq, double pamscale, FILE *fp, struct pstruct *ppst) {
+ int i, j, len;
+ int qi, rj;
+ int **pam2p;
+ int first, too_high;
+ char *query;
+ double freq, **freq2d, lambda, new_lambda;
+ double scale, scale_high, scale_low;
+
+ pam2p = ppst->pam2p[0];
+
+ if(1 != fread(&len, sizeof(int), 1, fp)) {
+ fprintf(stderr, "error reading from checkpoint file: %d\n", len);
+ exit(1);
+ }
+
+ if(len != n0) {
+ fprintf(stderr, "profile length (%d) and query length (%d) don't match!\n",
+ len,n0);
+ exit(1);
+ }
+
+ /* read over query sequence stored in BLAST profile */
+ if(NULL == (query = (char *) calloc(len, sizeof(char)))) {
+ fprintf(stderr, "Couldn't allocate memory for query!\n");
+ exit(1);
+ }
+
+ if(len != fread(query, sizeof(char), len, fp)) {
+ fprintf(stderr, "Couldn't read query sequence from profile: %s\n", query);
+ exit(1);
+ }
+
+ printf("%d\n%s\n",len,query);
+
+ /* currently we don't do anything with query; ideally, we should
+ check to see that it actually matches aa0 ... */
+
+ /* quick 2d array alloc: */
+ if((freq2d = (double **) calloc(n0, sizeof(double *))) == NULL) {
+ fprintf(stderr, "Couldn't allocate memory for frequencies!\n");
+ exit(1);
+ }
+
+ if((freq2d[0] = (double *) calloc(n0 * N_EFFECT, sizeof(double))) == NULL) {
+ fprintf(stderr, "Couldn't allocate memory for frequencies!\n");
+ exit(1);
+ }
+
+ /* a little pointer arithmetic to fill out 2d array: */
+ for (qi = 1 ; qi < n0 ; qi++) {
+ freq2d[qi] = freq2d[0] + (N_EFFECT * qi);
+ }
+
+ for (qi = 0 ; qi < n0 ; qi++) {
+ printf("%c",query[qi]);
+ for (rj = 0 ; rj < N_EFFECT ; rj++) {
+ if(1 != fread(&freq, sizeof(double), 1, fp)) {
+ fprintf(stderr, "Error while reading frequencies!\n");
+ exit(1);
+ }
+ printf(" %8.7g",freq*10.0);
+
+ if (freq > 1e-12) {
+ freq = log(freq /((double) (rrcounts[rj+1])/(double) rrtotal));
+ freq /= pamscale; /* this gets us close to originial pam scores */
+ freq2d[qi][rj] = freq;
+ }
+ else {freq2d[qi][rj] = freq;}
+ }
+ printf("\n");
+ }
+
+
+ /* now figure out the right scale */
+ scale = 1.0;
+ lambda = get_lambda(ppst->pam2[0], 20, 20, "\0ARNDCQEGHILKMFPSTWYV");
+
+ /* should be near 1.0 because of our initial scaling by ppst->pamscale */
+ fprintf(stderr, "real_lambda: %g\n", lambda);
+
+ /* get initial high/low scale values: */
+ first = 1;
+ while (1) {
+ fill_pam(pam2p, n0, 20, freq2d, scale);
+ new_lambda = get_lambda(pam2p, n0, 20, query);
+
+ if (new_lambda > lambda) {
+ if (first) {
+ first = 0;
+ scale = scale_high = 1.0 + 0.05;
+ scale_low = 1.0;
+ too_high = 1;
+ } else {
+ if (!too_high) break;
+ scale = (scale_high += scale_high - 1.0);
+ }
+ } else if (new_lambda > 0) {
+ if (first) {
+ first = 0;
+ scale_high = 1.0;
+ scale = scale_low = 1.0 - 0.05;
+ too_high = 0;
+ } else {
+ if (too_high) break;
+ scale = (scale_low += scale_low - 1.0);
+ }
+ } else {
+ fprintf(stderr, "new_lambda (%g) <= 0; matrix has positive average score", new_lambda);
+ exit(1);
+ }
+ }
+
+ /* now do binary search between low and high */
+ for (i = 0 ; i < 10 ; i++) {
+ scale = 0.5 * (scale_high + scale_low);
+ fill_pam(pam2p, n0, 20, freq2d, scale);
+ new_lambda = get_lambda(pam2p, n0, 20, query);
+
+ if (new_lambda > lambda) scale_low = scale;
+ else scale_high = scale;
+ }
+
+ scale = 0.5 * (scale_high + scale_low);
+ fill_pam(pam2p, n0, 20, freq2d, scale);
+
+ fprintf(stderr, "final scale: %g\n", scale);
+
+ for (qi = 0 ; qi < n0 ; qi++) {
+ fprintf(stderr, "%4d %c: ", qi+1, query[qi]);
+ for (rj = 1 ; rj <= 20 ; rj++) {
+ fprintf(stderr, "%4d", pam2p[qi][rj]);
+ }
+ fprintf(stderr, "\n");
+ }
+
+ free(freq2d[0]);
+ free(freq2d);
+
+ free(query);
+}
+
+/*
+ * alloc_pam(): allocates memory for the 2D pam matrix as well
+ * as for the integer array used to transmit the pam matrix
+ */
+void
+alloc_pam (int d1, int d2, struct pstruct *ppst)
+{
+ int i, *d2p;
+
+ if ((ppst->pam2[0] = (int **) malloc (d1 * sizeof (int *))) == NULL) {
+ fprintf(stderr,"Cannot allocate 2D pam matrix: %d",d1);
+ exit(1);
+ }
+
+ if ((ppst->pam2[1] = (int **) malloc (d1 * sizeof (int *))) == NULL) {
+ fprintf(stderr,"Cannot allocate 2D pam matrix: %d",d1);
+ exit(1);
+ }
+
+ if ((d2p = pam12 = (int *) malloc (d1 * d2 * sizeof (int))) == NULL) {
+ fprintf(stderr,"Cannot allocate 2D pam matrix: %d",d1);
+ exit(1);
+ }
+
+ for (i = 0; i < d1; i++, d2p += d2)
+ ppst->pam2[0][i] = d2p;
+
+ if ((d2p=pam12x= (int *) malloc (d1 * d2 * sizeof (int))) == NULL) {
+ fprintf(stderr,"Cannot allocate 2d pam matrix: %d",d2);
+ exit(1);
+ }
+
+ for (i = 0; i < d1; i++, d2p += d2)
+ ppst->pam2[1][i] = d2p;
+}
+
+void
+fill_pam(int **pam2p, int n0, int nsq, double **freq2d, double scale) {
+ int i, j;
+ double freq;
+
+ /* fprintf(stderr, "scale: %g\n", scale); */
+
+ /* now fill in the pam matrix: */
+ for (i = 0 ; i < n0 ; i++) {
+ for (j = 1 ; j <=nsq ; j++) {
+ freq = scale * freq2d[i][j-1];
+ if ( freq < 0.0) freq -= 0.5;
+ else freq += 0.5;
+ pam2p[i][j] = (int)(freq);
+ }
+ }
+}
+
+/*
+ * initpam2(struct pstruct pst): Converts 1-D pam matrix to 2-D
+ */
+void initpam2 (struct pstruct *ppst)
+{
+ int i, j, k, nsq, pam_xx, pam_xm;
+ int sa_x, sa_t, tmp;
+
+ nsq = ppst->nsq;
+ sa_x = pascii['X'];
+ sa_t = pascii['*'];
+
+ ppst->pam2[0][0][0] = -BIGNUM;
+ ppst->pam_h = -1; ppst->pam_l = 1;
+
+ k = 0;
+ for (i = 1; i <= nsq; i++) {
+ ppst->pam2[0][0][i] = ppst->pam2[0][i][0] = -BIGNUM;
+ for (j = 1; j <= i; j++) {
+ ppst->pam2[0][j][i] = ppst->pam2[0][i][j] = pam[k++] - ppst->pamoff;
+ if (ppst->pam_l > ppst->pam2[0][i][j]) ppst->pam_l =ppst->pam2[0][i][j];
+ if (ppst->pam_h < ppst->pam2[0][i][j]) ppst->pam_h =ppst->pam2[0][i][j];
+ }
+ }
+
+ ppst->nt_align = (ppst->dnaseq== SEQT_DNA || ppst->dnaseq == SEQT_RNA);
+
+ if (ppst->dnaseq == SEQT_RNA) {
+ tmp = ppst->pam2[0][nascii['G']][nascii['G']] - 1;
+ ppst->pam2[0][nascii['A']][nascii['G']] =
+ ppst->pam2[0][nascii['C']][nascii['T']] =
+ ppst->pam2[0][nascii['C']][nascii['U']] = tmp;
+ }
+
+ if (ppst->pam_x_set) {
+ for (i=1; i<=nsq; i++) {
+ ppst->pam2[0][sa_x][i] = ppst->pam2[0][i][sa_x]=ppst->pam_xm;
+ ppst->pam2[0][sa_t][i] = ppst->pam2[0][i][sa_t]=ppst->pam_xm;
+ }
+ ppst->pam2[0][sa_x][sa_x]=ppst->pam_xx;
+ ppst->pam2[0][sa_t][sa_t]=ppst->pam_xm;
+ }
+ else {
+ ppst->pam_xx = ppst->pam2[0][sa_x][sa_x];
+ ppst->pam_xm = ppst->pam2[0][1][sa_x];
+ }
+}
+
+double
+get_lambda(int **pam2p, int n0, int nsq, char *aa0) {
+ double lambda, H;
+ double *pr, tot, sum;
+ int i, ioff, j, min, max;
+
+ /* get min and max scores */
+ min = BIGNUM;
+ max = -BIGNUM;
+ if(pam2p[0][1] == -BIGNUM) {
+ ioff = 1;
+ n0++;
+ } else {
+ ioff = 0;
+ }
+
+ for (i = ioff ; i < n0 ; i++) {
+ for (j = 1; j <= nsq ; j++) {
+ if (min > pam2p[i][j])
+ min = pam2p[i][j];
+ if (max < pam2p[i][j])
+ max = pam2p[i][j];
+ }
+ }
+
+ /* fprintf(stderr, "min: %d\tmax:%d\n", min, max); */
+
+ if ((pr = (double *) calloc(max - min + 1, sizeof(double))) == NULL) {
+ fprintf(stderr, "Couldn't allocate memory for score probabilities: %d\n", max - min + 1);
+ exit(1);
+ }
+
+ tot = (double) rrtotal * (double) rrtotal * (double) n0;
+ for (i = ioff ; i < n0 ; i++) {
+ for (j = 1; j <= nsq ; j++) {
+ pr[pam2p[i][j] - min] +=
+ (double) ((double) rrcounts[aascii[aa0[i]]] * (double) rrcounts[j]) / tot;
+ }
+ }
+
+ sum = 0.0;
+ for(i = 0 ; i <= max-min ; i++) {
+ sum += pr[i];
+ /* fprintf(stderr, "%3d: %g %g\n", i+min, pr[i], sum); */
+ }
+ /* fprintf(stderr, "sum: %g\n", sum); */
+
+ for(i = 0 ; i <= max-min ; i++) { pr[i] /= sum; }
+
+ if (!karlin(min, max, pr, &lambda, &H)) {
+ fprintf(stderr, "Karlin lambda estimation failed\n");
+ }
+
+ /* fprintf(stderr, "lambda: %g\n", lambda); */
+ free(pr);
+
+ return lambda;
+}
+
+int **
+alloc_pam2p(int len, int nsq) {
+ int i;
+ int **pam2p;
+
+ if ((pam2p = (int **)calloc(len,sizeof(int *)))==NULL) {
+ fprintf(stderr," Cannot allocate pam2p: %d\n",len);
+ return NULL;
+ }
+
+ if((pam2p[0] = (int *)calloc((nsq+1)*len,sizeof(int)))==NULL) {
+ fprintf(stderr, "Cannot allocate pam2p[0]: %d\n", (nsq+1)*len);
+ free(pam2p);
+ return NULL;
+ }
+
+ for (i=1; i<len; i++) {
+ pam2p[i] = pam2p[0] + (i*(nsq+1));
+ }
+
+ return pam2p;
+}
+
+void free_pam2p(int **pam2p) {
+ if (pam2p) {
+ free(pam2p[0]);
+ free(pam2p);
+ }
+}
+
git://source.jalview.org / jabaws.git / blobdiff