forester/archive/RIO/others/hmmer/src/alphabet.c

   1 /************************************************************
   2  * HMMER - Biological sequence analysis with profile HMMs
   3  * Copyright (C) 1992-1999 Washington University School of Medicine
   4  * All Rights Reserved
   5  *
   6  *     This source code is distributed under the terms of the
   7  *     GNU General Public License. See the files COPYING and LICENSE
   8  *     for details.
   9  ************************************************************/
  10
  11 /* alphabet.c
  12  * Configuration of the global symbol alphabet information.
  13  * RCS $Id: alphabet.c,v 1.1.1.1 2005/03/22 08:34:08 cmzmasek Exp $
  14  */
  15
  16 #include <stdlib.h>
  17 #include <string.h>
  18 #include <ctype.h>
  19 #ifdef HMMER_THREADS
  20 #include <pthread.h>
  21 #endif /* HMMER_THREADS */
  22
  23 #include "config.h"
  24 #include "structs.h"
  25 #include "funcs.h"
  26 #include "squid.h"
  27
  28 #ifdef MEMDEBUG
  29 #include "dbmalloc.h"
  30 #endif
  31
  32 static void set_degenerate(char iupac, char *syms);
  33
  34
  35 /* Function: DetermineAlphabet()
  36  *
  37  * Purpose:  From a set of sequences (raw or aligned), make a good
  38  *           guess whether they're Nucleic, Amino, or something
  39  *           else, and set alphabet accordingly.
  40  *
  41  *           If Alphabet_type is already set, that means our
  42  *           autodetection was overridden from the command line,
  43  *           and we just set the other globals accordingly.
  44  */
  45 void
  46 DetermineAlphabet(char **rseqs, int  nseq)
  47 {
  48   int idx;
  49   int other, nucleic, amino;
  50   int type;
  51
  52   /* Autodetection of alphabet type.
  53    */
  54   type = hmmNOTSETYET;
  55   other = nucleic = amino = 0;
  56   for (idx = 0; idx < nseq; idx++) {
  57     switch (Seqtype(rseqs[idx])) {
  58     case kRNA:      nucleic++;   break;
  59     case kDNA:      nucleic++;   break;
  60     case kAmino:    amino++;     break;
  61     case kOtherSeq: other++;     break;
  62     default: Die("No such alphabet type");
  63     }
  64   }
  65
  66   if      (nucleic == nseq) type = hmmNUCLEIC;
  67   else if (amino   == nseq) type = hmmAMINO;
  68   else if (nucleic > amino && nucleic > other) {
  69     Warn("Looks like nucleic acid sequence, hope that's right");
  70     type = hmmNUCLEIC;
  71   }
  72   else if (amino > nucleic && amino > other) {
  73     Warn("Looks like amino acid sequence, hope that's right");
  74     type = hmmAMINO;
  75   }
  76   else Die("Sorry, I can't tell if that's protein or DNA");
  77
  78   /* Now set up the alphabet.
  79    */
  80   SetAlphabet(type);
  81 }
  82
  83
  84 /* Function: SetAlphabet()
  85  *
  86  * Purpose:  Set the alphabet globals, given an alphabet type
  87  *           of either hmmAMINO or hmmNUCLEIC.
  88  */
  89 void
  90 SetAlphabet(int type)
  91 {
  92   int x;
  93 #ifdef HMMER_THREADS
  94   pthread_mutex_t  alphabet_lock; /* alphabet is global; must protect to be threadsafe */
  95   int              rtn;           /* return code from pthreads */
  96
  97   if ((rtn = pthread_mutex_init(&alphabet_lock, NULL)) != 0)
  98     Die("pthread_mutex_init FAILED; %s\n", strerror(rtn));
  99   if ((rtn = pthread_mutex_lock(&alphabet_lock)) != 0)
 100     Die("pthread_mutex_lock FAILED: %s\n", strerror(rtn));
 101 #endif
 102
 103  /* Because the alphabet information is global, we must
 104   * be careful to make this a thread-safe function. The mutex
 105   * (above) takes care of that. But, indeed, it's also
 106   * just good sense (and more efficient) to simply never
 107   * allow resetting the alphabet. If type is Alphabet_type,
 108   * silently return; else die with an alphabet mismatch
 109   * warning.
 110   */
 111   if (Alphabet_type != hmmNOTSETYET)
 112     {
 113       if (type != Alphabet_type)
 114         Die("An alphabet type conflict occurred.\nYou probably mixed a DNA seq file with a protein model, or vice versa.");
 115
 116 #ifdef HMMER_THREADS
 117       if ((rtn = pthread_mutex_unlock(&alphabet_lock)) != 0)
 118         Die("pthread_mutex_unlock failure: %s\n", strerror(rtn));
 119 #endif
 120       return;
 121     }
 122
 123   switch(type) { /* Alphabet is not a string - careful! */
 124   case hmmAMINO:
 125     Alphabet_type     = type;
 126     strncpy(Alphabet, "ACDEFGHIKLMNPQRSTVWYBZX", 23);
 127     Alphabet_size     = 20;
 128     Alphabet_iupac    = 23;
 129     for (x = 0; x < Alphabet_iupac; x++) {
 130       memset(Degenerate[x], 0, Alphabet_size);
 131     }
 132     for (x = 0; x < Alphabet_size; x++) {
 133       Degenerate[x][x] = 1;
 134       DegenCount[x] = 1;
 135     }
 136     set_degenerate('B', "ND");
 137     set_degenerate('Z', "QE");
 138     set_degenerate('X', "ACDEFGHIKLMNPQRSTVWY");
 139     break;
 140   case hmmNUCLEIC:
 141     Alphabet_type     = type;
 142     strncpy(Alphabet, "ACGTUNRYMKSWHBVDX", 17);
 143     Alphabet_size     = 4;
 144     Alphabet_iupac    = 17;
 145     for (x = 0; x < Alphabet_iupac; x++) {
 146       memset(Degenerate[x], 0, Alphabet_size);
 147     }
 148     for (x = 0; x < Alphabet_size; x++) {
 149       Degenerate[x][x] = 1;
 150       DegenCount[x] = 1;
 151     }
 152     set_degenerate('U', "T");
 153     set_degenerate('N', "ACGT");
 154     set_degenerate('X', "ACGT");
 155     set_degenerate('R', "AG");
 156     set_degenerate('Y', "CT");
 157     set_degenerate('M', "AC");
 158     set_degenerate('K', "GT");
 159     set_degenerate('S', "CG");
 160     set_degenerate('W', "AT");
 161     set_degenerate('H', "ACT");
 162     set_degenerate('B', "CGT");
 163     set_degenerate('V', "ACG");
 164     set_degenerate('D', "AGT");
 165     break;
 166   default: Die("No support for non-nucleic or protein alphabets");
 167   }
 168
 169 #ifdef HMMER_THREADS
 170   if ((rtn = pthread_mutex_unlock(&alphabet_lock)) != 0)
 171     Die("pthread_mutex_unlock failure: %s\n", strerror(rtn));
 172 #endif
 173 }
 174
 175 /* Function: SymbolIndex()
 176  *
 177  * Purpose:  Convert a symbol to its index in Alphabet[].
 178  *           Bogus characters are converted to 'X'.
 179  *           More robust than the SYMIDX() macro but
 180  *           presumably slower.
 181  */
 182 int
 183 SymbolIndex(char sym)
 184 {
 185   char *s;
 186   return ((s = strchr(Alphabet, (char) toupper((int) sym))) == NULL) ?
 187           Alphabet_iupac-1 : s - Alphabet;
 188 }
 189
 190
 191 /* Function: DigitizeSequence()
 192  *
 193  * Purpose:  Internal representation of a sequence in HMMER is
 194  *           as a char array. 1..L are the indices
 195  *           of seq symbols in Alphabet[]. 0,L+1 are sentinel
 196  *           bytes, set to be Alphabet_iupac -- i.e. one more
 197  *           than the maximum allowed index.
 198  *
 199  *           Assumes that 'X', the fully degenerate character,
 200  *           is the last character in the allowed alphabet.
 201  *
 202  * Args:     seq - sequence to be digitized (0..L-1)
 203  *           L   - length of sequence
 204  *
 205  * Return:   digitized sequence, dsq.
 206  *           dsq is allocated here and must be free'd by caller.
 207  */
 208 char *
 209 DigitizeSequence(char *seq, int L)
 210 {
 211   char *dsq;
 212   int i;
 213
 214   dsq = MallocOrDie (sizeof(char) * (L+2));
 215   dsq[0] = dsq[L+1] = (char) Alphabet_iupac;
 216   for (i = 1; i <= L; i++)
 217     dsq[i] = SymbolIndex(seq[i-1]);
 218   return dsq;
 219 }
 220
 221
 222 /* Function: DedigitizeSequence()
 223  * Date:     SRE, Tue Dec 16 10:39:19 1997 [StL]
 224  *
 225  * Purpose:  Returns a 0..L-1 character string, converting the
 226  *           dsq back to the real alphabet.
 227  */
 228 char *
 229 DedigitizeSequence(char *dsq, int L)
 230 {
 231   char *seq;
 232   int i;
 233
 234   seq = MallocOrDie(sizeof(char) * (L+1));
 235   for (i = 0; i < L; i++)
 236     seq[i] = Alphabet[(int) dsq[i+1]];
 237   seq[L] = '\0';
 238   return seq;
 239 }
 240
 241
 242 /* Function: DigitizeAlignment()
 243  *
 244  * Purpose:  Given an alignment, return digitized unaligned
 245  *           sequence array. (Tracebacks are always relative
 246  *           to digitized unaligned seqs, even if they are
 247  *           faked from an existing alignment in modelmakers.c.)
 248  *
 249  * Args:     msa      - alignment to digitize
 250  *           ret_dsqs - RETURN: array of digitized unaligned sequences
 251  *
 252  * Return:   (void)
 253  *           dsqs is alloced here. Free2DArray(dseqs, nseq).
 254  */
 255 void
 256 DigitizeAlignment(MSA *msa, char ***ret_dsqs)
 257 {
 258   char **dsq;
 259   int    idx;                   /* counter for sequences     */
 260   int    dpos;                  /* position in digitized seq */
 261   int    apos;                  /* position in aligned seq   */
 262
 263   dsq = (char **) MallocOrDie (sizeof(char *) * msa->nseq);
 264   for (idx = 0; idx < msa->nseq; idx++) {
 265     dsq[idx] = (char *) MallocOrDie (sizeof(char) * (msa->alen+2));
 266
 267     dsq[idx][0] = (char) Alphabet_iupac; /* sentinel byte at start */
 268
 269     for (apos = 0, dpos = 1; apos < msa->alen; apos++) {
 270       if (! isgap(msa->aseq[idx][apos]))  /* skip gaps */
 271         dsq[idx][dpos++] = SymbolIndex(msa->aseq[idx][apos]);
 272     }
 273     dsq[idx][dpos] = (char) Alphabet_iupac; /* sentinel byte at end */
 274   }
 275   *ret_dsqs = dsq;
 276 }
 277
 278
 279 /* Function: P7CountSymbol()
 280  *
 281  * Purpose:  Given a possibly degenerate symbol code, increment
 282  *           a symbol counter array (generally an emission
 283  *           probability vector in counts form) appropriately.
 284  *
 285  * Args:     counters:  vector to count into. [0..Alphabet_size-1]
 286  *           symidx:    symbol index to count: [0..Alphabet_iupac-1]
 287  *           wt:        weight to use for the count; often 1.0
 288  *
 289  * Return:   (void)
 290  */
 291 void
 292 P7CountSymbol(float *counters, char symidx, float wt)
 293 {
 294   int x;
 295
 296   if (symidx < Alphabet_size)
 297     counters[(int) symidx] += wt;
 298   else
 299     for (x = 0; x < Alphabet_size; x++) {
 300       if (Degenerate[(int) symidx][x])
 301         counters[x] += wt / (float) DegenCount[(int) symidx];
 302     }
 303 }
 304
 305
 306 /* Function: DefaultGeneticCode()
 307  *
 308  * Purpose:  Configure aacode, mapping triplets to amino acids.
 309  *           Triplet index: AAA = 0, AAC = 1, ... UUU = 63.
 310  *           AA index: alphabetical: A=0,C=1... Y=19
 311  *           Stop codon: -1.
 312  *           Uses the stdcode1[] global translation table from SQUID.
 313  *
 314  * Args:     aacode  - preallocated 0.63 array for genetic code
 315  *
 316  * Return:   (void)
 317  */
 318 void
 319 DefaultGeneticCode(int *aacode)
 320 {
 321   int x;
 322
 323   for (x = 0; x < 64; x++) {
 324     if (*(stdcode1[x]) == '*') aacode[x] = -1;
 325     else                       aacode[x] = SYMIDX(*(stdcode1[x]));
 326   }
 327 }
 328
 329
 330 /* Function: DefaultCodonBias()
 331  *
 332  * Purpose:  Configure a codonbias table, mapping triplets to
 333  *           probability of using the triplet for the amino acid
 334  *           it represents: P(triplet | aa).
 335  *           The default is to assume codons are used equiprobably.
 336  *
 337  * Args:     codebias:  0..63 array of P(triplet|aa), preallocated.
 338  *
 339  * Return:   (void)
 340  */
 341 void
 342 DefaultCodonBias(float *codebias)
 343 {
 344   codebias[0]  = 1./2.; /* AAA Lys 2 */
 345   codebias[1]  = 1./2.; /* AAC Asn 2 */
 346   codebias[2]  = 1./2.; /* AAG Lys 2 */
 347   codebias[3]  = 1./2.; /* AAU Asn 2 */
 348   codebias[4]  = 1./4.; /* ACA Thr 4 */
 349   codebias[5]  = 1./4.; /* ACC Thr 4 */
 350   codebias[6]  = 1./4.; /* ACG Thr 4 */
 351   codebias[7]  = 1./4.; /* ACU Thr 4 */
 352   codebias[8]  = 1./6.; /* AGA Ser 6 */
 353   codebias[9]  = 1./6.; /* AGC Arg 6 */
 354   codebias[10] = 1./6.; /* AGG Ser 6 */
 355   codebias[11] = 1./6.; /* AGU Arg 6 */
 356   codebias[12] = 1./3.; /* AUA Ile 3 */
 357   codebias[13] = 1./3.; /* AUC Ile 3 */
 358   codebias[14] = 1.;    /* AUG Met 1 */
 359   codebias[15] = 1./3.; /* AUU Ile 3 */
 360   codebias[16] = 1./2.; /* CAA Gln 2 */
 361   codebias[17] = 1./2.; /* CAC His 2 */
 362   codebias[18] = 1./2.; /* CAG Gln 2 */
 363   codebias[19] = 1./2.; /* CAU His 2 */
 364   codebias[20] = 1./4.; /* CCA Pro 4 */
 365   codebias[21] = 1./4.; /* CCC Pro 4 */
 366   codebias[22] = 1./4.; /* CCG Pro 4 */
 367   codebias[23] = 1./4.; /* CCU Pro 4 */
 368   codebias[24] = 1./6.; /* CGA Arg 6 */
 369   codebias[25] = 1./6.; /* CGC Arg 6 */
 370   codebias[26] = 1./6.; /* CGG Arg 6 */
 371   codebias[27] = 1./6.; /* CGU Arg 6 */
 372   codebias[28] = 1./6.; /* CUA Leu 6 */
 373   codebias[29] = 1./6.; /* CUC Leu 6 */
 374   codebias[30] = 1./6.; /* CUG Leu 6 */
 375   codebias[31] = 1./6.; /* CUU Leu 6 */
 376   codebias[32] = 1./2.; /* GAA Glu 2 */
 377   codebias[33] = 1./2.; /* GAC Asp 2 */
 378   codebias[34] = 1./2.; /* GAG Glu 2 */
 379   codebias[35] = 1./2.; /* GAU Asp 2 */
 380   codebias[36] = 1./4.; /* GCA Ala 4 */
 381   codebias[37] = 1./4.; /* GCC Ala 4 */
 382   codebias[38] = 1./4.; /* GCG Ala 4 */
 383   codebias[39] = 1./4.; /* GCU Ala 4 */
 384   codebias[40] = 1./4.; /* GGA Gly 4 */
 385   codebias[41] = 1./4.; /* GGC Gly 4 */
 386   codebias[42] = 1./4.; /* GGG Gly 4 */
 387   codebias[43] = 1./4.; /* GGU Gly 4 */
 388   codebias[44] = 1./4.; /* GUA Val 4 */
 389   codebias[45] = 1./4.; /* GUC Val 4 */
 390   codebias[46] = 1./4.; /* GUG Val 4 */
 391   codebias[47] = 1./4.; /* GUU Val 4 */
 392   codebias[48] = 0.;    /* UAA och - */
 393   codebias[49] = 1./2.; /* UAC Tyr 2 */
 394   codebias[50] = 0.;    /* UAG amb - */
 395   codebias[51] = 1./2.; /* UAU Tyr 2 */
 396   codebias[52] = 1./6.; /* UCA Ser 6 */
 397   codebias[53] = 1./6.; /* UCC Ser 6 */
 398   codebias[54] = 1./6.; /* UCG Ser 6 */
 399   codebias[55] = 1./6.; /* UCU Ser 6 */
 400   codebias[56] = 0.;    /* UGA opa - */
 401   codebias[57] = 1./2.; /* UGC Cys 2 */
 402   codebias[58] = 1.;    /* UGG Trp 1 */
 403   codebias[59] = 1./2.; /* UGU Cys 2 */
 404   codebias[60] = 1./6.; /* UUA Leu 6 */
 405   codebias[61] = 1./2.; /* UUC Phe 2 */
 406   codebias[62] = 1./6.; /* UUG Leu 6 */
 407   codebias[63] = 1./2.; /* UUU Phe 2 */
 408 }
 409
 410
 411
 412 /* Function: set_degenerate()
 413  *
 414  * Purpose:  convenience function for setting up
 415  *           Degenerate[][] global for the alphabet.
 416  */
 417 static void
 418 set_degenerate(char iupac, char *syms)
 419 {
 420   DegenCount[strchr(Alphabet,iupac)-Alphabet] = strlen(syms);
 421   while (*syms) {
 422     Degenerate[strchr(Alphabet,iupac)-Alphabet]
 423               [strchr(Alphabet,*syms)-Alphabet] = 1;
 424     syms++;
 425   }
 426 }