binaries/src/fasta34/pssm_asn_subs.c

   1 /* pssm_asn_subs.c */
   2
   3
   4 /* $Name: fa_34_26_5 $ - $Id: pssm_asn_subs.c,v 1.15 2007/04/02 18:08:11 wrp Exp $ */
   5
   6 /* copyright (C) 2005 by William R. Pearson and the U. of Virginia */
   7
   8 /* this code is designed to parse the ASN.1 binary encoded scoremat
   9    object produced by blastpgp -C file.ckpt_asn -u 2 */
  10
  11 #include <stdio.h>
  12 #include <stdlib.h>
  13 #include <string.h>
  14
  15 #include "defs.h"
  16
  17 int parse_pssm_asn();
  18 int parse_pssm2_asn();
  19
  20 int
  21 parse_pssm_asn_fa(FILE *afd, int *n_rows, int *n_cols,
  22                   unsigned char **query, double ***freqs,
  23                   char *matrix, int *gap_open, int *gap_extend,
  24                   double *lambda);
  25
  26
  27
  28 #define COMPO_NUM_TRUE_AA 20
  29
  30 /**positions of true characters in protein alphabet*/
  31 /*
  32 static int trueCharPositions[COMPO_NUM_TRUE_AA] = {
  33   1,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,22
  34 };
  35 */
  36
  37 #define COMPO_LARGEST_ALPHABET 28
  38
  39 /*
  40 static char ncbieaatoa[COMPO_LARGEST_ALPHABET] = {"-ABCDEFGHIJKLMNOPQRSTUVWXYZ"};
  41
  42 static int alphaConvert[COMPO_LARGEST_ALPHABET] = {
  43   (-1), 0, (-1), 4, 3, 6, 13, 7, 8, 9, 11, 10, 12, 2, 14, 5, 1, 15,
  44   16, 19,   17, (-1), 18, (-1), (-1), (-1), (-1), (-1)
  45 };
  46 */
  47
  48 int pssm_aa_order[20] = { 1,  /*A*/
  49                           16, /*R*/
  50                           13, /*N*/
  51                            4, /*D*/
  52                            3, /*C*/
  53                           15, /*Q*/
  54                            5, /*E*/
  55                            7, /*G*/
  56                            8, /*H*/
  57                            9, /*I*/
  58                           11, /*L*/
  59                           10, /*K*/
  60                           12, /*M*/
  61                            6, /*F*/
  62                           14, /*P*/
  63                           17, /*S*/
  64                           18, /*T*/
  65                           20, /*W*/
  66                           22, /*Y*/
  67                           19}; /*V*/
  68
  69
  70 #define ASN_SEQ 48
  71 #define ASN_SEQOF 49
  72
  73 #define ASN_PSSM_QUERY 166
  74 #define ASN_PSSM2_QUERY 162
  75
  76 #define ASN_PSSM_IS_PROT 160
  77 #define ASN_PSSM2_MATRIX 161
  78 #define ASN_PSSM_NROWS 162
  79 #define ASN_PSSM_NCOLS 163
  80
  81 #define ASN_PSSM2_NCOLS 163
  82 #define ASN_PSSM2_NROWS 164
  83 #define ASN_PSSM_BYCOL 165
  84 #define ASN_PSSM_INTERMED_DATA 167
  85 #define ASN_PSSM_FREQS 162
  86 #define ASN_PSSM2_FREQS 165
  87 #define ASN_PSSM2_LAMBDA 166
  88
  89 #define ASN_IS_STR 26
  90 #define ASN_IS_INT  2
  91 #define ASN_IS_BOOL 1
  92 #define ASN_IS_ENUM 10
  93
  94 struct asn_bstruct {
  95   FILE *fd;
  96   unsigned char *buf;
  97   unsigned char *abp;
  98   unsigned char *buf_max;
  99   int len;
 100 };
 101
 102 #define ASN_BUF 1024
 103
 104 unsigned char *
 105 chk_asn_buf(struct asn_bstruct *asnp, int v) {
 106   int new_buf;
 107
 108   if (v > ASN_BUF) {
 109     fprintf(stderr," attempt to read %d bytes ASN.1 data > buffer size (%d)\n",
 110             v, ASN_BUF);
 111     exit(1);
 112   }
 113
 114   if (asnp->abp + v > asnp->buf_max) {
 115
 116     /* move down the left over stuff */
 117     asnp->len = asnp->buf_max - asnp->abp;
 118
 119     memmove(asnp->buf, asnp->abp, asnp->len);
 120
 121     asnp->abp = asnp->buf;
 122     new_buf = ASN_BUF - asnp->len;
 123
 124     if (!feof(asnp->fd) &&
 125         (new_buf=fread(asnp->buf + asnp->len, sizeof(char), new_buf, asnp->fd)) != 0) {
 126       asnp->len += new_buf;
 127     }
 128
 129     asnp->buf_max = asnp->buf + asnp->len;
 130
 131     if (asnp->len < v) {
 132       fprintf(stderr, " Unable to read %d bytes\n",v);
 133       exit(1);
 134     }
 135   }
 136   /* otherwise, v bytes are currently in the buffer */
 137
 138   return asnp->abp;
 139 }
 140
 141 /* read_asn_dest reads v bytes into oct_str if v <= o_len */
 142 /* read_asn_dest is required for ASN data entities that are longer than ASN_BUF (1024) */
 143 unsigned char *
 144 read_asn_dest(struct asn_bstruct *asnp, int v, unsigned char *oct_str, int o_len) {
 145   int new_buf;
 146   unsigned char *oct_ptr;
 147
 148
 149   if (v > o_len) {
 150     fprintf(stderr, " read_asn_dest - cannot read %d bytes into %d buffer\n",
 151             v, o_len);
 152     exit(1);
 153   }
 154
 155   if (asnp->abp + v <= asnp->buf_max) {
 156     memmove(oct_str, asnp->abp, v);
 157     return asnp->abp+v;
 158   }
 159   else {
 160     /* move down the left over stuff */
 161
 162     asnp->len = asnp->buf_max - asnp->abp;
 163
 164     memmove(oct_str, asnp->abp, asnp->len);
 165     oct_ptr = oct_str+asnp->len;
 166     v -= asnp->len;
 167
 168     asnp->abp = asnp->buf;
 169     new_buf = ASN_BUF;
 170
 171     while ((new_buf=fread(asnp->buf, sizeof(char), new_buf, asnp->fd)) != 0) {
 172       asnp->len = new_buf;
 173       asnp->buf_max = asnp->buf + asnp->len;
 174       if (v <= new_buf) {       /* we have it all this time */
 175         memmove(oct_ptr, asnp->buf, v);
 176         asnp->len -= v;
 177         asnp->abp = asnp->buf + v;
 178         break;
 179       }
 180       else {    /* we need to read some more */
 181         memmove(oct_ptr, asnp->buf, new_buf);
 182         v -= new_buf;
 183         new_buf = ASN_BUF;
 184       }
 185     }
 186   }
 187   return asnp->buf + v;
 188 }
 189
 190 unsigned char *
 191 get_astr_bool(struct asn_bstruct *asnp, int *val) {
 192
 193   int v_len, v;
 194
 195   asnp->abp = chk_asn_buf(asnp,5);
 196
 197   v = 0;
 198   if (*asnp->abp++ != 1) { /* check for int */
 199     fprintf(stderr," bool missing\n");
 200   }
 201   else {
 202     v_len = *asnp->abp++;
 203     if (v_len != 1) {
 204       fprintf(stderr, "boolean length != 1 : %d\n", v_len);
 205       v = *asnp->abp++;
 206     }
 207     else { v = *asnp->abp++;}
 208   }
 209   asnp->abp += 2;       /* skip over null's */
 210   *val = v;
 211   return asnp->abp;
 212 }
 213
 214 unsigned char *
 215 get_astr_int(struct asn_bstruct *asnp,
 216             int *val) {
 217
 218   int v_len, v;
 219
 220   v = 0;
 221
 222   asnp->abp = chk_asn_buf(asnp,8);
 223
 224   if (*asnp->abp++ != 2) { /* check for int */
 225     fprintf(stderr," int missing\n");
 226   }
 227   else {
 228     v_len = *asnp->abp++;
 229     while (v_len-- > 0) {
 230       v *= 256;
 231       v += *asnp->abp++;
 232     }
 233     asnp->abp += 2;     /* skip over null's */
 234   }
 235   *val = v;
 236   return asnp->abp;
 237 }
 238
 239 unsigned char *
 240 get_astr_enum(struct asn_bstruct *asnp, int *val) {
 241
 242   int v_len, v;
 243
 244   asnp->abp = chk_asn_buf(asnp,5);
 245
 246   v = 0;
 247   if (*asnp->abp++ != ASN_IS_ENUM) { /* check for int */
 248     fprintf(stderr," enum missing\n");
 249   }
 250   else {
 251     v_len = *asnp->abp++;
 252     while (v_len-- > 0) { v *= 256;  v += *asnp->abp++; }
 253     asnp->abp += 2;     /* skip over null's */
 254   }
 255   *val = v;
 256
 257   return asnp->abp;
 258 }
 259
 260 unsigned char *
 261 get_astr_packedfloat(struct asn_bstruct *asnp, double *val) {
 262
 263   int v_len, v;
 264   char tmp_str[64];
 265
 266   asnp->abp = chk_asn_buf(asnp,2);
 267
 268   v = 0;
 269   if (*asnp->abp++ != 9) { /* check for packed float */
 270     fprintf(stderr," float missing\n");
 271     *val = 0;
 272     return asnp->abp;
 273   }
 274   else {
 275     v_len = *asnp->abp++;
 276
 277     if (v_len > 63) {
 278       fprintf(stderr," real string too long: %d\n",v_len);
 279     }
 280
 281     asnp->abp = chk_asn_buf(asnp,v_len);
 282
 283     if (v_len == 2  && *asnp->abp == '\0' && *(asnp->abp+1)=='0') {
 284       asnp->abp += 2;
 285       *val = 0.0;
 286     }
 287     else {      /* copy and scan it */
 288       if (*asnp->abp != '\0') {
 289         fprintf(stderr, " packedfloat - expected 0, got %d\n", *asnp->abp);
 290         *val = -1.0;
 291         return asnp->abp;
 292       }
 293       asnp->abp++;
 294       strncpy(tmp_str, (char *)asnp->abp, sizeof(tmp_str)-1);
 295       tmp_str[v_len-1] = '\0';
 296       tmp_str[63] = '\0';
 297       sscanf(tmp_str,"%lg",val);
 298       asnp->abp += v_len-1;
 299     }
 300   }
 301   return asnp->abp;
 302 }
 303
 304 unsigned char *
 305 get_astr_str(struct asn_bstruct *asnp, char *text, int t_len) {
 306
 307   int v_len;
 308
 309   asnp->abp = chk_asn_buf(asnp,2);
 310
 311   text[0] = '\0';
 312   if (*asnp->abp++ != ASN_IS_STR) { /* check for str */
 313     fprintf(stderr," str missing\n");
 314   }
 315   else {
 316     v_len = *asnp->abp++;
 317     if (v_len > 128) { /* need to read the length from the next bytes */
 318       t_len = v_len &0x7f;
 319
 320       asnp->abp = chk_asn_buf(asnp,t_len);
 321
 322       for (v_len =0; t_len; t_len--) { v_len = (v_len << 8) + *asnp->abp++; }
 323     }
 324
 325     /* read v_len bytes */
 326
 327     asnp->abp = read_asn_dest(asnp,v_len, (unsigned char *)text, t_len);
 328     asnp->abp += 2; /* skip over last nulls */
 329   }
 330   return asnp->abp;
 331 }
 332
 333 #define ASN_BIOSEQ_SEQ 160
 334 #define ASN_BIOSEQ_ID  160
 335 #define ASN_BIOSEQ_ID_VAL 160
 336
 337 #define ASN_BIOSEQ_ID_LOCAL 161
 338 #define ASN_BIOSEQ_ID_GIBBSQ 162
 339 #define ASN_BIOSEQ_ID_GIBBMT 163
 340 #define ASN_BIOSEQ_ID_GB 164
 341 #define ASN_BIOSEQ_ID_EMBL 165
 342 #define ASN_BIOSEQ_ID_PIR 166
 343 #define ASN_BIOSEQ_ID_SP 167
 344 #define ASN_BIOSEQ_ID_PATENT 168
 345 #define ASN_BIOSEQ_ID_OTHER 169
 346 #define ASN_BIOSEQ_ID_GEN 170
 347 #define ASN_BIOSEQ_ID_GI 171
 348
 349 #define ASN_BIOSEQ_TEXTID_NAME 160
 350 #define ASN_BIOSEQ_TEXTID_ACC 161
 351 #define ASN_BIOSEQ_TEXTID_REL 162
 352 #define ASN_BIOSEQ_TEXTID_VER 163
 353
 354 #define ASN_BIOSEQ_DESCR  161
 355 #define ASN_BIOSEQ_INST  162
 356 #define ASN_BIOSEQ_TITLE  164
 357 #define ASN_BIOSEQ_INST_REPR  160
 358 #define ASN_BIOSEQ_INST_MOL  161
 359 #define ASN_BIOSEQ_INST_LEN  162
 360 #define ASN_BIOSEQ_INST_TOPOL  166
 361 #define ASN_BIOSEQ_INST_SEQD  167
 362 #define ASN_OCTET_STR  65
 363 #define ASN_NCBIeaa 65
 364
 365 unsigned char *
 366 get_astr_seqdescr(struct asn_bstruct *asnp,
 367                  char *descr) {
 368
 369   int end_seq=0;
 370
 371   /* get seqof '1' */
 372   /* get 164/128 -  title */
 373   /* get string */
 374   /* pop nulls */
 375
 376   asnp->abp = chk_asn_buf(asnp,6);
 377
 378   if (*asnp->abp == ASN_SEQOF) {
 379     end_seq++;
 380     asnp->abp += 2;
 381   }
 382   else {
 383     fprintf(stderr, " missing ASN_SEQOF '1': %0x %0x\n",*asnp->abp, asnp->abp[1]);
 384   }
 385
 386   if (*asnp->abp == ASN_BIOSEQ_TITLE) {
 387     asnp->abp+=2;
 388     asnp->abp = get_astr_str(asnp, descr, MAX_STR);
 389   }
 390   else {
 391     fprintf(stderr, " missing ASN_BIOSEQ_TITLE '1': %0x %0x\n",*asnp->abp, asnp->abp[1]);
 392   }
 393
 394   asnp->abp = chk_asn_buf(asnp,2);
 395
 396   asnp->abp += 2; /* skip over nulls */
 397
 398   return asnp->abp;
 399 }
 400
 401 unsigned char *
 402 get_astr_octstr(struct asn_bstruct *asnp,
 403                unsigned char *oct_str,
 404                int o_len) {
 405
 406   int q_len, v_len;
 407
 408   asnp->abp = chk_asn_buf(asnp,2);
 409
 410   if (*asnp->abp++ == ASN_NCBIeaa) {
 411     /* get length  of length */
 412     if (*asnp->abp > 128) {
 413       v_len = *asnp->abp++ & 0x7f;
 414
 415       asnp->abp = chk_asn_buf(asnp,v_len);
 416
 417       q_len = 0;
 418       while (v_len-- > 0) {
 419         q_len *= 256;
 420         q_len += *asnp->abp++;
 421       }
 422     }
 423     else {
 424       q_len = *asnp->abp++ & 0x7f;
 425     }
 426
 427     asnp->abp = read_asn_dest(asnp, q_len, oct_str, o_len);
 428
 429     oct_str[min(q_len,o_len)]='\0';
 430
 431     asnp->abp += 2;     /* skip characters and NULL's */
 432   }
 433   return asnp->abp;
 434 }
 435
 436 unsigned char *
 437 get_astr_seqinst(struct asn_bstruct *asnp,
 438                 unsigned char **query,
 439                 int *nq) {
 440
 441   int end_seq=0, tmp;
 442
 443   /* get sequence '0' */
 444   /* get 160/128/10/len/val -  repr enum raw val */
 445   /* get 161/128/10/len/val -  mol enum aa val */
 446   /* get 162/128/02/len/val -  length int val */
 447   /* get 166/128 - topology (empty) */
 448   /* get 167/128 - seq-data */
 449   /* get 65/len+128/len/octet_string */
 450   /* pop nulls */
 451
 452   asnp->abp = chk_asn_buf(asnp,12);
 453
 454   if (*asnp->abp == ASN_SEQ) {
 455     end_seq++;
 456     asnp->abp += 2;
 457   }
 458   else {
 459     fprintf(stderr, " missing ASN_SEQ '0': %0x %0x\n",*asnp->abp, asnp->abp[1]);
 460   }
 461
 462   if (*asnp->abp == ASN_BIOSEQ_INST_REPR && *(asnp->abp+1) == 128) {
 463     asnp->abp+=2;
 464     asnp->abp = get_astr_enum(asnp, &tmp);
 465   }
 466   else {
 467     fprintf(stderr, " missing ASN_BIOSEQ_INST_REPR 160: %0x %0x\n",*asnp->abp, asnp->abp[1]);
 468   }
 469
 470   if (*asnp->abp == ASN_BIOSEQ_INST_MOL && *(asnp->abp+1) == 128) {
 471     asnp->abp+=2;
 472     asnp->abp = get_astr_enum(asnp, &tmp);
 473   }
 474   else {
 475     fprintf(stderr, " missing ASN_BIOSEQ_INST_MOL 161: %0x %0x\n",*asnp->abp, asnp->abp[1]);
 476   }
 477
 478   if (*asnp->abp == ASN_BIOSEQ_INST_LEN) {
 479     asnp->abp+=2;
 480     asnp->abp = get_astr_int(asnp, nq);
 481   }
 482   else {
 483     fprintf(stderr, " missing ASN_BIOSEQ_INST_LEN 161: %0x %0x\n",*asnp->abp, asnp->abp[1]);
 484     return asnp->abp;
 485   }
 486
 487   if ((*query = (unsigned char *)calloc(*nq + 1, sizeof(char)))==NULL) {
 488     fprintf(stderr, " cannot read %d char query\n", *nq+1);
 489   }
 490
 491   if (*asnp->abp == ASN_BIOSEQ_INST_TOPOL && *(asnp->abp+1) == 128 ) {
 492     asnp->abp += 2;
 493   }
 494
 495   if (*asnp->abp == ASN_BIOSEQ_INST_SEQD) {
 496     asnp->abp+=2;
 497     asnp->abp = get_astr_octstr(asnp, *query, *nq );
 498   }
 499   else {
 500     fprintf(stderr, " missing ASN_BIOSEQ_INST_SEQD 166: %0x %0x\n",*asnp->abp, asnp->abp[1]);
 501     return asnp->abp;
 502   }
 503
 504   asnp->abp += 4; /* skip over nulls */
 505
 506   return asnp->abp;
 507 }
 508
 509
 510 unsigned char *
 511 get_astr_textid( struct asn_bstruct *asnp,
 512                 char *name,
 513                 char *acc) {
 514   int end_seq = 0;
 515   int ver;
 516
 517   chk_asn_buf(asnp,16);
 518
 519   if (*asnp->abp != ASN_SEQ) {
 520     fprintf(stderr, " Expected ASN_SEQ: %0x %0x\n",*asnp->abp, asnp->abp[1]);
 521   }
 522   else {asnp->abp += 2; end_seq++;}
 523
 524   name[0] = acc[0] = '\0';
 525
 526   while (*asnp->abp != '\0') {
 527     if (*asnp->abp == ASN_BIOSEQ_TEXTID_NAME) {
 528       asnp->abp+=2;
 529       asnp->abp = get_astr_str(asnp, name, MAX_SSTR);
 530     }
 531     if (*asnp->abp == ASN_BIOSEQ_TEXTID_ACC) {
 532       asnp->abp+=2;
 533       asnp->abp = get_astr_str(asnp, acc, MAX_SSTR);
 534     }
 535     if (*asnp->abp == ASN_BIOSEQ_TEXTID_VER) {
 536       asnp->abp+=2;
 537       asnp->abp = get_astr_int(asnp, &ver);
 538     }
 539   }
 540   asnp->abp += 4;
 541   while (end_seq-- > 0) { asnp->abp += 4; }
 542   return asnp->abp;
 543 }
 544
 545 unsigned char *
 546 get_astr_query(struct asn_bstruct *asnp,
 547               int *gi,
 548               char *name,
 549               char *acc,
 550               char *descr,
 551               unsigned char **query,
 552               int *nq
 553               ) {
 554
 555   int end_seq = 0;
 556
 557   asnp->abp = chk_asn_buf(asnp,32);
 558
 559   if (*asnp->abp != ASN_BIOSEQ_SEQ) {
 560     fprintf(stderr, "Bioseq - missing SEQ 1: %2x %2x\n",*asnp->abp, asnp->abp[1]);
 561     return asnp->abp;
 562   }
 563   else { asnp->abp += 2;}
 564
 565   if (*asnp->abp != ASN_SEQ && *asnp->abp != ASN_SEQOF ) {
 566     fprintf(stderr, "Bioseq - missing SEQUENCE tag 1: %2x %2x\n",*asnp->abp, asnp->abp[1]);
 567     return asnp->abp;
 568   }
 569   else {
 570     end_seq++;
 571     asnp->abp += 2;
 572   }
 573
 574   if (*asnp->abp != ASN_BIOSEQ_ID) {
 575     fprintf(stderr, "Bioseq - missing ID tag: %2x %2x\n",*asnp->abp, asnp->abp[1]);
 576     return asnp->abp;
 577   }
 578   else {
 579     asnp->abp += 2;
 580     if (*asnp->abp != ASN_SEQOF) {
 581       fprintf(stderr, "missing bioseq/id/SEQOF tag: %d\n",*asnp->abp);
 582       return asnp->abp;
 583     }
 584     else {
 585       asnp->abp += 2;
 586       if (*asnp->abp == ASN_BIOSEQ_ID_VAL && *(asnp->abp+1)==128) { asnp->abp += 2;}
 587
 588       if (*asnp->abp == ASN_BIOSEQ_ID_GI ) {
 589         asnp->abp+=2;
 590         asnp->abp = get_astr_int(asnp, gi);
 591       }
 592
 593       if (*asnp->abp == ASN_BIOSEQ_ID_LOCAL) {
 594         *gi = 0;
 595         acc[0] = '\0';
 596
 597         asnp->abp+=2;
 598         asnp->abp = get_astr_str(asnp, name, MAX_SSTR);
 599         asnp->abp += 2;
 600       }
 601       else if (*asnp->abp == ASN_BIOSEQ_ID_SP || *asnp->abp == ASN_BIOSEQ_ID_EMBL ||
 602           *asnp->abp == ASN_BIOSEQ_ID_GB || *asnp->abp == ASN_BIOSEQ_ID_PIR ||
 603           *asnp->abp == ASN_BIOSEQ_ID_OTHER ) {
 604
 605         asnp->abp+=2;
 606         asnp->abp = get_astr_textid(asnp, name, acc);
 607       }
 608     }
 609   }
 610
 611   while (*asnp->abp == 0) asnp->abp += 2;
 612
 613   if (*asnp->abp == ASN_BIOSEQ_DESCR) {
 614     asnp->abp+=2;
 615     asnp->abp = get_astr_seqdescr(asnp, descr);
 616     asnp->abp += 2;             /* skip nulls */
 617   }
 618   else { descr[0] = '\0';}
 619
 620   if (*asnp->abp != ASN_BIOSEQ_INST) {
 621     fprintf(stderr, "Bioseq - missing ID tag: %2x %2x\n",*asnp->abp, asnp->abp[1]);
 622     return asnp->abp;
 623   }
 624   else {
 625     asnp->abp += 2;
 626     asnp->abp = get_astr_seqinst(asnp, query, nq);
 627     asnp->abp += 2;             /* skip nulls */
 628   }
 629   return asnp->abp;
 630 }
 631
 632 unsigned char *
 633 get_astr_query2(struct asn_bstruct *asnp,
 634               int *gi,
 635               char *name,
 636               char *acc,
 637               char *descr,
 638               unsigned char **query,
 639               int *nq
 640               ) {
 641
 642   int end_seq = 0;
 643
 644   asnp->abp = chk_asn_buf(asnp,32);
 645
 646   if (*asnp->abp != ASN_BIOSEQ_SEQ) {
 647     fprintf(stderr, "Bioseq - missing SEQ 1: %2x %2x\n",*asnp->abp, asnp->abp[1]);
 648     return asnp->abp;
 649   }
 650   else { asnp->abp += 2;}
 651
 652   if (*asnp->abp != ASN_SEQOF ) {
 653     fprintf(stderr, "Bioseq2 - missing SEQOF tag 1: %2x %2x\n",*asnp->abp, asnp->abp[1]);
 654     return asnp->abp;
 655   }
 656   else {
 657     end_seq++;
 658     asnp->abp += 2;
 659   }
 660
 661   if (*asnp->abp != ASN_BIOSEQ_ID) {
 662     fprintf(stderr, "Bioseq - missing ID tag: %2x %2x\n",*asnp->abp, asnp->abp[1]);
 663     return asnp->abp;
 664   }
 665   else {
 666     asnp->abp += 2;
 667     if (*asnp->abp == ASN_SEQOF) {
 668       asnp->abp += 2;
 669     }
 670
 671     if (*asnp->abp == ASN_BIOSEQ_ID_VAL && *(asnp->abp+1)==128) { asnp->abp += 2;}
 672
 673     if (*asnp->abp == ASN_BIOSEQ_ID_GI ) {
 674       asnp->abp+=2;
 675       asnp->abp = get_astr_int(asnp, gi);
 676     }
 677
 678     if (*asnp->abp == ASN_BIOSEQ_ID_LOCAL) {
 679       *gi = 0;
 680       acc[0] = '\0';
 681
 682       asnp->abp+=2;
 683       asnp->abp = get_astr_str(asnp, name, MAX_SSTR);
 684       asnp->abp += 2;
 685       }
 686     else if (*asnp->abp == ASN_BIOSEQ_ID_SP || *asnp->abp == ASN_BIOSEQ_ID_EMBL ||
 687              *asnp->abp == ASN_BIOSEQ_ID_GB || *asnp->abp == ASN_BIOSEQ_ID_PIR ||
 688              *asnp->abp == ASN_BIOSEQ_ID_OTHER ) {
 689
 690       asnp->abp+=2;
 691       asnp->abp = get_astr_textid(asnp, name, acc);
 692     }
 693   }
 694
 695   while (*asnp->abp == 0) asnp->abp += 2;
 696
 697   if (*asnp->abp == ASN_BIOSEQ_DESCR) {
 698     asnp->abp+=2;
 699     asnp->abp = get_astr_seqdescr(asnp, descr);
 700     asnp->abp += 2;             /* skip nulls */
 701   }
 702   else { descr[0] = '\0';}
 703
 704   if (*asnp->abp != ASN_BIOSEQ_INST) {
 705     fprintf(stderr, "Bioseq - missing ID tag: %2x %2x\n",*asnp->abp, asnp->abp[1]);
 706     return asnp->abp;
 707   }
 708   else {
 709     asnp->abp += 2;
 710     asnp->abp = get_astr_seqinst(asnp, query, nq);
 711     asnp->abp += 2;             /* skip nulls */
 712   }
 713   return asnp->abp;
 714 }
 715
 716 unsigned char *
 717 get_pssm_freqs(struct asn_bstruct *asnp,
 718                double **freqs,
 719                int n_rows,
 720                int n_cols,
 721                int by_row) {
 722
 723   int i_rows, i_cols;
 724   int in_seq = 0;
 725
 726   double f_val;
 727
 728   asnp->abp = chk_asn_buf(asnp,4);
 729
 730   if (*asnp->abp == ASN_SEQ) {
 731     in_seq = 1;
 732     asnp->abp += 2;
 733     in_seq = 1;
 734   }
 735
 736   if (!by_row) {
 737     for (i_cols = 0; i_cols < n_cols; i_cols++) {
 738       for (i_rows = 0; i_rows < n_rows; i_rows++) {
 739         asnp->abp = get_astr_packedfloat(asnp, &f_val);
 740         freqs[i_cols][i_rows] = f_val;
 741       }
 742     }
 743   }
 744   else {
 745     for (i_rows = 0; i_rows < n_rows; i_rows++) {
 746       for (i_cols = 0; i_cols < n_cols; i_cols++) {
 747         asnp->abp = get_astr_packedfloat(asnp, &f_val);
 748         freqs[i_cols][i_rows] = f_val;
 749       }
 750     }
 751   }
 752   if (in_seq) {asnp->abp +=2;}  /* skip nulls */
 753   asnp->abp += 2;
 754   return asnp->abp;
 755 }
 756
 757 unsigned char *
 758 get_pssm_intermed(struct asn_bstruct *asnp,
 759                   double **freqs,
 760                   int n_rows,
 761                   int n_cols,
 762                   int by_row) {
 763
 764   asnp->abp = chk_asn_buf(asnp,4);
 765
 766   if (*asnp->abp == ASN_SEQ) {
 767     asnp->abp += 2;
 768     if (*asnp->abp == ASN_PSSM_FREQS) {
 769       asnp->abp+=2;
 770       asnp->abp = get_pssm_freqs(asnp, freqs, n_rows, n_cols, by_row);
 771     }
 772     asnp->abp +=2;      /* skip nulls */
 773   }
 774   asnp->abp += 2;
 775   return asnp->abp;
 776 }
 777
 778
 779 #define ASN_PSSM_PARAMS 161
 780 #define ASN_PSSM_PARAMS_PSEUDOCNT 160
 781 #define ASN_PSSM_PARAMS_RPSPARAMS 161
 782 #define ASN_PSSM_RPSPARAMS_MATRIX 160
 783 #define ASN_PSSM_RPSPARAMS_GAPOPEN 161
 784 #define ASN_PSSM_RPSPARAMS_GAPEXT 162
 785
 786 unsigned char *
 787 get_pssm_rpsparams(struct asn_bstruct *asnp,
 788                char *matrix,
 789                int *gap_open,
 790                int *gap_ext) {
 791
 792   int end_seq=0;
 793
 794   asnp->abp = chk_asn_buf(asnp,4);
 795
 796   if (*asnp->abp == ASN_SEQ) {
 797     asnp->abp += 2;
 798     end_seq++;
 799   }
 800
 801   if (*asnp->abp == ASN_PSSM_RPSPARAMS_MATRIX) {
 802     asnp->abp+=2;
 803     asnp->abp = get_astr_str(asnp, matrix, MAX_SSTR);
 804   }
 805
 806   if (*asnp->abp == ASN_PSSM_RPSPARAMS_GAPOPEN) {
 807     asnp->abp+=2;
 808     asnp->abp = get_astr_int(asnp, gap_open);
 809   }
 810
 811   if (*asnp->abp == ASN_PSSM_RPSPARAMS_GAPEXT) {
 812     asnp->abp+=2;
 813     asnp->abp = get_astr_int(asnp, gap_ext);
 814   }
 815
 816   if (end_seq) { chk_asn_buf(asnp,end_seq * 2); }
 817   while (end_seq-- > 0) { asnp->abp += 2; }
 818   return asnp->abp;
 819 }
 820
 821 unsigned char *
 822 get_pssm_params(struct asn_bstruct *asnp,
 823                int *pseudo_cnts,
 824                char *matrix,
 825                int *gap_open,
 826                int *gap_ext) {
 827
 828   int end_seq=0;
 829
 830   asnp->abp = chk_asn_buf(asnp,6);
 831
 832   if (*asnp->abp == ASN_SEQ) {
 833     asnp->abp += 2;
 834     end_seq++;
 835   }
 836
 837   if (*asnp->abp == ASN_PSSM_PARAMS_PSEUDOCNT) {
 838     asnp->abp+=2;
 839     asnp->abp = get_astr_int(asnp, pseudo_cnts);
 840   }
 841
 842   if (*asnp->abp == ASN_PSSM_PARAMS_RPSPARAMS) {
 843     asnp->abp+=2;
 844     asnp->abp = get_pssm_rpsparams(asnp, matrix, gap_open, gap_ext);
 845     asnp->abp += 2;
 846   }
 847   while (end_seq-- > 0) { asnp->abp+=2; }
 848   return asnp->abp;
 849 }
 850
 851
 852 unsigned char *
 853 get_pssm2_intermed(struct asn_bstruct *asnp,
 854                    double ***freqs,
 855                    int n_rows,
 856                    int n_cols) {
 857
 858   int i;
 859   double **my_freqs;
 860
 861   if ((my_freqs = (double **) calloc(n_cols, sizeof(double *)))==NULL) {
 862     fprintf(stderr, " cannot allocate freq cols - %d\n", n_cols);
 863     exit(1);
 864   }
 865
 866   if ((my_freqs[0] = (double *) calloc(n_cols * n_rows, sizeof(double)))==NULL) {
 867     fprintf(stderr, " cannot allocate freq rows * cols - %d * %d\n", n_rows, n_cols);
 868     exit(1);
 869   }
 870
 871   for (i=1; i < n_cols; i++) {
 872     my_freqs[i] = my_freqs[i-1] + n_rows;
 873   }
 874
 875   *freqs = my_freqs;
 876
 877   chk_asn_buf(asnp, 8);
 878
 879   return get_pssm_freqs(asnp, my_freqs, n_rows, n_cols, 0);
 880 }
 881
 882 int
 883 parse_pssm2_asn(struct asn_bstruct *asnp,
 884                 int *gi,
 885                 char *name,
 886                 char *acc,
 887                 char *descr,
 888                 unsigned char **query,
 889                 int *nq,
 890                 int *n_rows,
 891                 int *n_cols,
 892                 double ***freqs,
 893                 int *pseudo_cnts,
 894                 char *matrix,
 895                 double *lambda_p) {
 896
 897   int is_protein;
 898   int have_rows, have_cols;
 899
 900   chk_asn_buf(asnp, 32);
 901
 902   if (memcmp(asnp->abp, "\241\2000\200",4) != 0) {
 903     fprintf(stderr, "improper PSSM2 start\n");
 904     return -1;
 905   }
 906   else {asnp->abp+=4;}
 907
 908   if (*asnp->abp == ASN_BIOSEQ_SEQ ) {
 909     asnp->abp = get_astr_query2(asnp, gi, name, acc, descr, query, nq);
 910   }
 911
 912   /* finish up the nulls */
 913   while (*asnp->abp == '\0') { asnp->abp += 2;}
 914
 915   if (*asnp->abp == ASN_PSSM2_QUERY &&
 916         asnp->abp[2] != ASN_SEQ ) {
 917       fprintf(stderr, "improper PSSM2 start\n");
 918       return -1;
 919   }
 920   else {asnp->abp += 4;}
 921
 922   while (*asnp->abp != '\0' ) {
 923
 924     switch (*asnp->abp) {
 925     case ASN_PSSM_IS_PROT :
 926       asnp->abp+=2;
 927       asnp->abp = get_astr_bool(asnp, &is_protein);
 928       break;
 929
 930     case ASN_PSSM2_MATRIX :
 931       asnp->abp+=2;
 932       asnp->abp = get_astr_str(asnp, matrix, MAX_SSTR);
 933       break;
 934
 935     case ASN_PSSM2_NROWS :
 936       asnp->abp+=2;
 937       asnp->abp = get_astr_int(asnp, n_rows);
 938
 939       if (*n_rows > 0) { have_rows = 1; }
 940       else {
 941         fprintf(stderr, " bad n_row count\n");
 942         exit(1);
 943       }
 944       break;
 945
 946     case  ASN_PSSM2_NCOLS :
 947       asnp->abp+=2;
 948       asnp->abp = get_astr_int(asnp, n_cols);
 949       if (*n_cols > 0) {
 950         have_cols = 1;
 951       }
 952       else {
 953         fprintf(stderr, " bad n_row count\n");
 954         exit(1);
 955       }
 956       break;
 957
 958     case ASN_PSSM2_FREQS :
 959       asnp->abp += 4;
 960       if (*asnp->abp == '\0') { asnp->abp += 4;}
 961       break;
 962
 963     case ASN_PSSM2_LAMBDA :
 964       asnp->abp += 2;
 965       asnp->abp = get_astr_packedfloat(asnp,lambda_p);
 966       asnp->abp +=2;    /* skip over end of ASN_PSSM2_LAMBDA */
 967       break;
 968
 969     case ASN_PSSM_INTERMED_DATA :
 970       asnp->abp += 2;
 971       asnp->abp = get_pssm2_intermed(asnp, freqs, *n_rows, *n_cols);
 972       asnp->abp += 4;
 973       break;
 974
 975     default: asnp->abp += 2;
 976     }
 977   }
 978
 979
 980   return 1;
 981 }
 982
 983 int
 984 parse_pssm_asn(FILE *afd,
 985                int *gi,
 986                char *name,
 987                char *acc,
 988                char *descr,
 989                unsigned char **query,
 990                int *nq,
 991                int *n_rows,
 992                int *n_cols,
 993                double ***freqs,
 994                int *pseudo_cnts,
 995                char *matrix,
 996                int *gap_open,
 997                int *gap_ext,
 998                double *lambda_p) {
 999
1000   int is_protein, pssm_version;
1001   int i;
1002   int have_rows, have_cols, by_col;
1003   double **my_freqs;
1004
1005   struct asn_bstruct asn_str;
1006
1007   if ((asn_str.buf = (unsigned char *)calloc(ASN_BUF, sizeof(char))) == NULL ) {
1008     fprintf(stderr, " cannot allocate asn_buf (%d)\n",ASN_BUF);
1009     exit(1);
1010   }
1011
1012   asn_str.fd = afd;
1013   asn_str.len = ASN_BUF;
1014   asn_str.abp = asn_str.buf_max = asn_str.buf + ASN_BUF;
1015
1016   chk_asn_buf(&asn_str, 32);
1017
1018   if (memcmp(asn_str.abp, "0\200\240\200",4) != 0) {
1019     fprintf(stderr, "improper PSSM header -");
1020     return -1;
1021   }
1022   else {asn_str.abp+=4;}
1023
1024   if (*asn_str.abp == ASN_IS_INT) {
1025     asn_str.abp = get_astr_int(&asn_str, &pssm_version);
1026     if (pssm_version != 2) {
1027       fprintf(stderr, "PSSM2 version mismatch: %d\n",pssm_version);
1028       return -1;
1029     }
1030     *gap_open = *gap_ext = 0;
1031     return parse_pssm2_asn(&asn_str, gi, name, acc, descr,
1032                            query, nq,
1033                            n_rows, n_cols, freqs,
1034                            pseudo_cnts, matrix,
1035                            lambda_p);
1036   }
1037
1038   if (*asn_str.abp == ASN_SEQ) { asn_str.abp += 2;  }
1039
1040   if (*asn_str.abp == ASN_PSSM_IS_PROT ) {
1041     asn_str.abp+=2;
1042     asn_str.abp = get_astr_bool(&asn_str, &is_protein);
1043   }
1044
1045   if (*asn_str.abp == ASN_PSSM_NROWS ) {
1046     asn_str.abp+=2;
1047     asn_str.abp = get_astr_int(&asn_str, n_rows);
1048
1049     if (*n_rows > 0) { have_rows = 1; }
1050     else {
1051       fprintf(stderr, " bad n_row count\n");
1052       exit(1);
1053     }
1054   }
1055
1056   if (*asn_str.abp == ASN_PSSM_NCOLS ) {
1057     asn_str.abp+=2;
1058     asn_str.abp = get_astr_int(&asn_str, n_cols);
1059     if (*n_cols > 0) {
1060       have_cols = 1;
1061     }
1062     else {
1063       fprintf(stderr, " bad n_row count\n");
1064       exit(1);
1065     }
1066   }
1067
1068   if (*asn_str.abp == ASN_PSSM_BYCOL ) {
1069     asn_str.abp+=2;
1070     asn_str.abp = get_astr_bool(&asn_str, &by_col);
1071   }
1072
1073   /* we have read everything up to the query
1074
1075      n_cols gives us the query length, which we can allocate;
1076   */
1077
1078   if (*asn_str.abp == ASN_PSSM_QUERY ) {
1079     asn_str.abp+=2;
1080     asn_str.abp = get_astr_query(&asn_str, gi, name, acc, descr, query, nq);
1081     *nq = *n_cols;
1082   }
1083
1084   /* finish up the nulls */
1085   while (*asn_str.abp == '\0') { asn_str.abp += 2;}
1086
1087   if (*asn_str.abp == ASN_PSSM_INTERMED_DATA) {
1088
1089     if (!have_rows || !have_cols) {
1090       fprintf(stderr, " cannot allocate freq - missing rows/cols - %d/%d\n",
1091               have_rows, have_cols);
1092       return -1;
1093     }
1094
1095     if ((my_freqs = (double **) calloc(*n_cols, sizeof(double *)))==NULL) {
1096       fprintf(stderr, " cannot allocate freq cols - %d\n", *n_cols);
1097       return -1;
1098     }
1099
1100     if ((my_freqs[0] = (double *) calloc(*n_cols * *n_rows, sizeof(double)))==NULL) {
1101       fprintf(stderr, " cannot allocate freq rows * cols - %d * %d\n", *n_rows, *n_cols);
1102       return -1;
1103     }
1104     for (i=1; i < *n_cols; i++) {
1105       my_freqs[i] = my_freqs[i-1] + *n_rows;
1106     }
1107
1108     *freqs = my_freqs;
1109
1110     asn_str.abp+=2;
1111     asn_str.abp = get_pssm_intermed(&asn_str, my_freqs, *n_rows, *n_cols, by_col);
1112     asn_str.abp += 4;
1113   }
1114
1115   if (*asn_str.abp == ASN_PSSM_PARAMS ) {
1116       asn_str.abp+=2;
1117       asn_str.abp = get_pssm_params(&asn_str, pseudo_cnts, matrix, gap_open, gap_ext);
1118   }
1119   else if (*asn_str.abp == 0) {asn_str.abp+=2;}
1120   return 1;
1121 }
1122
1123 int
1124 parse_pssm_asn_fa( FILE *fd,
1125                    int *n_rows_p, int *n_cols_p,
1126                    unsigned char **query,
1127                    double ***freq2d,
1128                    char *matrix,
1129                    int *gap_open_p,
1130                    int *gap_extend_p,
1131                    double *lambda_p) {
1132
1133   int qi, rj;
1134   int gi;
1135   double tmp_freqs[COMPO_LARGEST_ALPHABET];
1136   char name[MAX_SSTR], acc[MAX_SSTR], descr[MAX_STR];
1137   int  nq;
1138   int pseudo_cnts, ret_val;
1139
1140   /* parse the file */
1141
1142   ret_val = parse_pssm_asn(fd, &gi, name, acc, descr, query, &nq,
1143                            n_rows_p, n_cols_p, freq2d,
1144                            &pseudo_cnts, matrix, gap_open_p, gap_extend_p,
1145                            lambda_p);
1146
1147   if (ret_val <=0) return ret_val;
1148
1149   /* transform the frequencies */
1150
1151   for (qi = 0; qi < *n_cols_p; qi++) {
1152     for (rj = 0; rj < *n_rows_p; rj++) { tmp_freqs[rj] = (*freq2d)[qi][rj];}
1153
1154     for (rj = 0; rj < COMPO_NUM_TRUE_AA; rj++) {
1155       (*freq2d)[qi][rj] = tmp_freqs[pssm_aa_order[rj]];
1156     }
1157   }
1158   return 1;
1159 }