binaries/src/fasta34/apam.c

   1 /*      pam.c   19-June-86
   2         copyright (c) 1987 William R. Pearson
   3         read in the alphabet and pam matrix data
   4         designed for universal matcher
   5
   6         This version reads BLAST format (square) PAM files
   7 */
   8
   9 /* $Name: fa_34_26_5 $ - $Id: apam.c,v 1.41 2007/03/31 18:47:20 wrp Exp $ */
  10
  11 #include <stdio.h>
  12 #include <stdlib.h>
  13 #include <string.h>
  14 #include <ctype.h>
  15
  16 #include "defs.h"
  17 #include "param.h"
  18
  19 #define XTERNAL
  20 #include "uascii.h"
  21 #include "upam.h"
  22 #undef XTERNAL
  23
  24 extern void alloc_pam (int d1, int d2, struct pstruct *ppst);
  25
  26 void
  27 pam_opts(char *smstr, struct pstruct *ppst) {
  28   char *bp;
  29
  30   ppst->pam_ms = 0;
  31   ppst->pamoff = 0;
  32
  33   if ((bp=strchr(smstr,'-'))!=NULL) {
  34     if (!strncmp(bp+1,"MS",2) || !strncmp(bp+1,"ms",2)) {
  35       ppst->pam_ms = 1;
  36     }
  37     else {
  38       ppst->pamoff=atoi(bp+1);
  39     }
  40     *bp = '\0';
  41    }
  42   else if ((bp=strchr(smstr,'+'))!=NULL) {
  43     ppst->pamoff= -atoi(bp+1);
  44     *bp = '\0';
  45   }
  46 }
  47
  48 /* modified 13-Oct-2005 to accomodate assymetrical matrices */
  49
  50 int
  51 initpam (char *mfname, struct pstruct *ppst)
  52 {
  53    char    line[512], *lp;
  54    int     i, j, iaa, pval;
  55    int *hsq, nsq;
  56    int *sascii;
  57    char *sq;
  58    int ess_tmp, max_val, min_val;
  59    int have_es = 0;
  60    FILE   *fmat;
  61
  62    pam_opts(mfname, ppst);
  63
  64    if ((fmat = fopen (mfname, "r")) == NULL)
  65    {
  66       printf ("***WARNING*** cannot open scoring matrix file %s\n", mfname);
  67       fprintf (stderr,"***WARNING*** cannot open scoring matrix file %s\n", mfname);
  68       return 0;
  69    }
  70
  71 /*
  72    the size of the alphabet is determined in advance
  73 */
  74    hsq = ppst->hsq;
  75    sq = ppst->sq;
  76
  77    ppst->nt_align = (ppst->dnaseq == SEQT_DNA || ppst->dnaseq == SEQT_RNA);
  78
  79 /*
  80      look for alphabet line, skipping the comments
  81      alphabet ends up in line[]
  82 */
  83    while (fgets (line, sizeof(line), fmat) != NULL && line[0]=='#');
  84
  85    /* decide whether this is a protein or DNA matrix */
  86    if (ppst->nt_align) sascii = &nascii[0];
  87    else sascii = &aascii[0];
  88
  89 /*
  90   re-initialize sascii[] for matrix alphabet
  91 */
  92
  93    /* save ',' value used by FASTS/FASTM/FASTF */
  94    ess_tmp = sascii[','];
  95
  96 /*      clear out sascii        */
  97    for (i = 0; i <= AAMASK; i++) sascii[i] = NA;
  98
  99 /*      set end of line stop    */
 100    sascii[0] = sascii['\r'] = sascii['\n'] = EL;
 101
 102    sascii[','] = ess_tmp;
 103
 104 /* read the alphabet - determine alphabet nsq */
 105    sq[0] = '\0';
 106    for (i = 0, nsq = 1; line[i]; i++) {
 107      if (line[i] == '*') have_es = 1;
 108      if (line[i] > ' ') sq[nsq++] = toupper (line[i]);
 109    }
 110    sq[nsq]='\0';
 111    nsq--;
 112
 113 /* set end of sequence stop */
 114    fprintf(stderr,"sq[%d]: %s\n",nsq,sq+1);
 115
 116 /* initialize sascii */
 117    for (iaa = 1; iaa <= nsq; iaa++) {
 118      sascii[sq[iaa]] = iaa;
 119    }
 120    if (ppst->dnaseq==SEQT_DNA) {
 121      sascii['U'] = sascii['T'];
 122      sascii['u'] = sascii['t'];
 123    }
 124    else if (ppst->dnaseq==SEQT_RNA) {
 125      sascii['T'] = sascii['U'];
 126      sascii['t'] = sascii['u'];
 127    }
 128
 129 /*
 130    finished with sascii[]
 131 */
 132
 133 /*
 134   setup hnt (ambiguous nt hash) values
 135 */
 136    hsq[0] = 0;
 137    for (iaa = 1; iaa <= nsq; iaa++) {
 138      hsq[iaa]=iaa;
 139    }
 140    if (ppst->nt_align) {        /* DNA ambiguitities */
 141      hsq[sascii['R']]=hsq[sascii['M']]=hsq[sascii['W']]=hsq[sascii['A']];
 142      hsq[sascii['D']]=hsq[sascii['H']]=hsq[sascii['V']]=hsq[sascii['A']];
 143      hsq[sascii['N']]=hsq[sascii['X']]=hsq[sascii['A']];
 144      hsq[sascii['Y']]=hsq[sascii['S']]=hsq[sascii['B']]=hsq[sascii['C']];
 145      hsq[sascii['K']]=hsq[sascii['G']];
 146    }
 147    else         /* protein ambiguities */
 148      if (ppst->dnaseq == SEQT_UNK || ppst->dnaseq == SEQT_PROT ||
 149          (ppst->nsq >= 20 && ppst->nsq <= 24)) {
 150      hsq[sascii['B']] = hsq[sascii['N']];
 151      hsq[sascii['Z']] = hsq[sascii['E']];
 152      hsq[sascii['X']] = hsq[sascii['A']];
 153    }
 154    /* here if non-DNA, non-protein sequence */
 155    else ppst->dnaseq = SEQT_OTHER;
 156
 157 /*
 158   check for 2D pam  - if not found, allocate it
 159 */
 160
 161    if (!ppst->have_pam2) {
 162      alloc_pam (MAXSQ, MAXSQ, ppst);
 163      ppst->have_pam2 = 1;
 164    }
 165
 166 /*
 167   read the scoring matrix values
 168 */
 169
 170    max_val = -1;
 171    min_val =  1;
 172    for (j=0; j < nsq; j++) ppst->pam2[0][0][j] = -BIGNUM;
 173    for (iaa = 1; iaa <= nsq; iaa++) {   /* read pam value line */
 174      if (fgets(line,sizeof(line),fmat)==NULL) {
 175        fprintf (stderr," error reading pam line: %s\n",line);
 176        exit (1);
 177      }
 178      /*     fprintf(stderr,"%d/%d %s",iaa,nsq,line); */
 179      strtok(line," \t\n");              /* skip the letter (residue) */
 180      ppst->pam2[0][i][0] = -BIGNUM;
 181      for (j = 1; j <= nsq; j++) {       /* iaa limits to triangle */
 182        lp=strtok(NULL," \t\n");         /* get the number string */
 183        pval=ppst->pam2[0][iaa][j]=atoi(lp);     /* convert to integer */
 184        if (pval > max_val) max_val = pval;
 185        if (pval < min_val) min_val = pval;
 186      }
 187    }
 188
 189    if (have_es==0) {
 190      sascii['*']=nsq;
 191      nsq++;
 192      sq[nsq]='*';
 193      sq[nsq+1]='\0';
 194      for (j=1; j<=nsq; j++) ppst->pam2[0][nsq][j]= -1;
 195      ppst->pam2[0][nsq][nsq]= max_val/2;
 196    }
 197
 198    ppst->sqx[0]='\0';   /* initialize sqx[] */
 199    for (i=1; i<= nsq; i++) {
 200      ppst->sqx[i] = sq[i];
 201      ppst->sqx[i+nsq] = tolower(sq[i]);
 202      if (sascii[aa[i]] < NA && sq[i] >= 'A' && sq[i] <= 'Z')
 203        sascii[aa[i] - 'A' + 'a'] = sascii[aa[i]]+nsq;
 204    }
 205
 206    ppst->nsq = nsq;     /* save new nsq */
 207    ppst->nsqx = nsq*2;  /* save new nsqx */
 208
 209    ppst->pam_h = max_val;
 210    ppst->pam_l = min_val;
 211
 212    strncpy (ppst->pamfile, mfname, MAX_FN);
 213    ppst->pamfile[MAX_FN-1]='\0';
 214
 215    if (ppst->pam_ms) {
 216      strncat(ppst->pamfile,"-MS",MAX_FN-strlen(ppst->pamfile)-1);
 217    }
 218    ppst->pamfile[MAX_FN-1]='\0';
 219    fclose (fmat);
 220    return 1;
 221 }
 222
 223 /* make a DNA scoring from +match/-mismatch values */
 224
 225 void mk_n_pam(int *arr,int siz, int mat, int mis)
 226 {
 227   int i, j, k;
 228   /* current default match/mismatch values */
 229   int max_mat = +5;
 230   int min_mis = -4;
 231   float f_val, f_scale;
 232
 233   f_scale = (float)(mat - mis)/(float)(max_mat - min_mis);
 234
 235   k = 0;
 236   for (i = 0; i<nnt-1; i++)
 237     for (j = 0; j <= i; j++ ) {
 238       if (arr[k] == max_mat) arr[k] = mat;
 239       else if (arr[k] == min_mis) arr[k] = mis;
 240       else if (arr[k] != -1) {
 241         f_val = (arr[k] - min_mis)*f_scale + 0.5;
 242         arr[k] = f_val + mis;
 243       }
 244       k++;
 245     }
 246 }
 247
 248 struct std_pam_str {
 249   char abbrev[6];
 250   char name[10];
 251   int *pam;
 252   float scale;
 253   int gdel, ggap;
 254 };
 255
 256 static
 257 struct std_pam_str std_pams[] = {
 258   {"P120", "PAM120", apam120, 0.346574, -20, -3},
 259   {"P250", "PAM250", apam250, 0.231049, -12, -2},
 260   {"P10",  "MD10",  a_md10,  0.346574, -27, -4},
 261   {"M10",  "MD10",  a_md10,  0.346574, -27, -4},
 262   {"MD10", "MD10",  a_md10,  0.346574, -27, -4},
 263   {"P20",  "MD20",  a_md20,  0.346574, -26, -4},
 264   {"M20",  "MD20",  a_md20,  0.346574, -26, -4},
 265   {"MD20", "MD20",  a_md20,  0.346574, -26, -4},
 266   {"P40",  "MD40",  a_md40,  0.346574, -25, -4},
 267   {"M40",  "MD40",  a_md40,  0.346574, -25, -4},
 268   {"MD40", "MD40",  a_md40,  0.346574, -25, -4},
 269   {"BL50", "BL50",   abl50,  0.231049, -12, -2},
 270   {"BL62", "BL62",   abl62,  0.346574,  -8, -1},
 271   {"BP62", "BL62",   abl62,  0.346574, -12, -1},
 272   {"BL80", "BL80",   abl80,  0.346574, -12, -2},
 273   {"\0",   "\0",     NULL,   0.0,        0,  0}
 274 };
 275
 276 int
 277 standard_pam(char *smstr, struct pstruct *ppst, int del_set, int gap_set) {
 278
 279   struct std_pam_str *std_pam_p;
 280
 281   pam_opts(smstr, ppst);
 282
 283   for (std_pam_p = std_pams; std_pam_p->abbrev[0]; std_pam_p++ ) {
 284     if (strcmp(smstr,std_pam_p->abbrev)==0) {
 285       pam = std_pam_p->pam;
 286       strncpy(ppst->pamfile,std_pam_p->name,MAX_FN);
 287       ppst->pamfile[MAX_FN-1]='\0';
 288       if (ppst->pam_ms) {
 289         strncat(ppst->pamfile,"-MS",MAX_FN-strlen(ppst->pamfile)-1);
 290       }
 291       ppst->pamfile[MAX_FN-1]='\0';
 292 #ifdef OLD_FASTA_GAP
 293       if (!del_set) ppst->gdelval = std_pam_p->gdel;
 294 #else
 295       if (!del_set) ppst->gdelval = std_pam_p->gdel-std_pam_p->ggap;
 296 #endif
 297       if (!gap_set) ppst->ggapval = std_pam_p->ggap;
 298       ppst->pamscale = std_pam_p->scale;
 299       return 1;
 300     }
 301   }
 302   return 0;
 303 }
 304
 305 /* ESS must match uascii.h */
 306 #define ESS 49
 307
 308 void
 309 build_xascii(int *qascii, char *save_str) {
 310   int i, max_save;
 311   int comma_val, term_val;
 312   int save_arr[MAX_SSTR];
 313
 314   comma_val = qascii[','];
 315   term_val = qascii['*'];
 316
 317   /* preserve special characters */
 318   for (i=0; i < MAX_SSTR && save_str[i]; i++ ) {
 319     save_arr[i] = qascii[save_str[i]];
 320   }
 321   max_save = i;
 322
 323   for (i=1; i<128; i++) {
 324     qascii[i]=NA;
 325   }
 326   /* range of values in aax, ntx is from 1..naax,nntx -
 327      do not zero-out qascii[0] - 9 Oct 2002 */
 328
 329   for (i=1; i<naax; i++) {
 330     qascii[aax[i]]=aax[i];
 331   }
 332
 333   for (i=1; i<nntx; i++) {
 334     qascii[ntx[i]]=ntx[i];
 335   }
 336
 337   qascii['\n']=qascii['\r']=qascii[0] = EL;
 338
 339   qascii[','] = comma_val;
 340   qascii['*'] = term_val;
 341
 342   for (i=0; i < max_save; i++) {
 343     qascii[save_str[i]]=save_arr[i];
 344   }
 345 }
 346
 347 /*
 348    checks for lower case letters in *sq array;
 349    if not present, map lowercase to upper
 350 */
 351 void
 352 init_ascii(int is_ext, int *sascii, int is_dna) {
 353
 354   int isq, have_lc;
 355   char *sq, term_char;
 356   int nsq;
 357
 358   if (is_dna==SEQT_UNK) return;
 359
 360   term_char = sascii['*'];
 361
 362   if (is_dna==SEQT_DNA || is_dna == SEQT_RNA) {
 363     if (is_ext) {
 364       sq = &ntx[0];
 365       nsq = nntx;
 366     }
 367     else {sq = &nt[0]; nsq = nnt;}
 368   }
 369   else {
 370     if (is_ext) { sq = &aax[0]; nsq = naax; }
 371     else {sq = &aa[0]; nsq = naa;}
 372   }
 373
 374
 375 /* initialize sascii from sq[], checking for lower-case letters */
 376   have_lc = 0;
 377   for (isq = 1; isq <= nsq; isq++) {
 378      sascii[sq[isq]] = isq;
 379      if (sq[isq] >= 'a' && sq[isq] <= 'z') have_lc = 1;
 380   }
 381
 382   /* no lower case letters in alphabet, map lower case to upper */
 383   if (have_lc != 1) {
 384     for (isq = 1; isq <= nsq; isq++) {
 385       if (sq[isq] >= 'A' && sq[isq] <= 'Z') sascii[sq[isq]-'A'+'a'] = isq;
 386     }
 387     if (is_dna==1) sascii['u'] = sascii['t'];
 388   }
 389
 390   sascii['*']=term_char;
 391 }
 392
 393 print_pam(struct pstruct *ppst) {
 394   int i, nsq, ip;
 395   char *sq;
 396
 397   fprintf(stderr," ext_sq_set: %d\n",ppst->ext_sq_set);
 398
 399   nsq = ppst->nsq;
 400   ip = 0;
 401   sq = ppst->sq;
 402
 403   fprintf(stderr," sq[%d]: %s\n",nsq, sq);
 404
 405   if (ppst->ext_sq_set) {
 406     nsq = ppst->nsqx;
 407     ip = 1;
 408     sq = ppst->sqx;
 409     fprintf(stderr," sq[%d]: %s\n",nsq, sq);
 410   }
 411
 412   for (i=1; i<=nsq; i++) {
 413     fprintf(stderr," %c:%c - %3d\n",sq[i], sq[i], ppst->pam2[ip][i][i]);
 414   }
 415 }