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

   1 /************************************************************
   2  * HMMER - Biological sequence analysis with profile-HMMs
   3  * Copyright (C) 1992-1997 Sean R. Eddy
   4  *
   5  *   This source code is distributed under the terms of the
   6  *   GNU General Public License. See the files COPYING and
   7  *   GNULICENSE for details.
   8  *
   9  ************************************************************/
  10
  11 /* states.c
  12  *
  13  * alloc, free, and initialization of state structures
  14  */
  15
  16 #include <stdio.h>
  17 #include <stdlib.h>
  18 #include <string.h>
  19 #include <math.h>
  20 #include "squid.h"
  21 #include "config.h"
  22 #include "structs.h"
  23 #include "funcs.h"
  24
  25 #ifdef MEMDEBUG
  26 #include "dbmalloc.h"
  27 #endif
  28
  29
  30 struct hmm_struc *
  31 AllocHMM(int M)                         /* length of model to make */
  32 {
  33   struct hmm_struc *hmm;        /* RETURN: blank HMM */
  34
  35   hmm        = (struct hmm_struc *)   MallocOrDie (sizeof(struct hmm_struc));
  36   hmm->ins   = (struct basic_state *) MallocOrDie (sizeof(struct basic_state) * (M+2));
  37   hmm->del   = (struct basic_state *) MallocOrDie (sizeof(struct basic_state) * (M+2));
  38   hmm->mat   = (struct basic_state *) MallocOrDie (sizeof(struct basic_state) * (M+2));
  39   hmm->ref   = (char *)  MallocOrDie ((M+2) * sizeof(char));
  40   hmm->cs    = (char *)  MallocOrDie ((M+2) * sizeof(char));
  41   hmm->xray  = (float *) MallocOrDie ((M+2) * sizeof(float) * NINPUTS);
  42   hmm->M     = M;
  43   hmm->name  = Strdup("unnamed"); /* name is not optional. */
  44
  45   hmm->flags = 0;
  46   ZeroHMM(hmm);
  47   return hmm;
  48 }
  49
  50 /* Function: ZeroHMM()
  51  *
  52  * Purpose:  Zero emission and transition counts in an HMM.
  53  */
  54 void
  55 ZeroHMM(struct hmm_struc *hmm)
  56 {
  57   int k, ts, idx;
  58
  59   for (k = 0; k <= hmm->M+1; k++)
  60     {
  61       for (ts = 0; ts < 3; ts++)
  62         {
  63           hmm->mat[k].t[ts] = 0.0;
  64           hmm->ins[k].t[ts] = 0.0;
  65           hmm->del[k].t[ts] = 0.0;
  66         }
  67       for (idx = 0; idx < Alphabet_size; idx++)
  68         {
  69           hmm->mat[k].p[idx]   = 0.0;
  70           hmm->ins[k].p[idx]   = 0.0;
  71           hmm->del[k].p[idx]   = 0.0;
  72         }
  73     }
  74 }
  75
  76
  77 /* Function: LogifyHMM()
  78  *
  79  * Purpose:  Convert a probability-form HMM to log probabilities.
  80  *           Best to do this on a modifiable copy of an HMM.
  81  */
  82 void
  83 LogifyHMM(struct hmm_struc *hmm)
  84 {
  85   int k, ts, idx;
  86
  87   for (k = 0; k <= hmm->M+1; k++)
  88     {
  89       for (ts = 0; ts < 3; ts++)
  90         {
  91           hmm->mat[k].t[ts] = sreLOG2(hmm->mat[k].t[ts]);
  92           hmm->ins[k].t[ts] = sreLOG2(hmm->ins[k].t[ts]);
  93           hmm->del[k].t[ts] = sreLOG2(hmm->del[k].t[ts]);
  94         }
  95       for (idx = 0; idx < Alphabet_size; idx++)
  96         {
  97           hmm->mat[k].p[idx]   = sreLOG2(hmm->mat[k].p[idx]);
  98           hmm->ins[k].p[idx]   = sreLOG2(hmm->ins[k].p[idx]);
  99         }
 100     }
 101 }
 102
 103 /* Function: LogoddsifyHMM()
 104  *
 105  * Convert a probability form HMM to log odds scores.
 106  * Best to do this on a modifiable copy of an HMM.
 107  */
 108 void
 109 LogoddsifyHMM(struct hmm_struc *hmm)
 110 {
 111   int k, ts, x;
 112
 113   for (k = 0; k <= hmm->M+1; k++)
 114     {
 115       for (ts = 0; ts < 3; ts++)
 116         {
 117           hmm->mat[k].t[ts] = sreLOG2(hmm->mat[k].t[ts]);
 118           hmm->ins[k].t[ts] = sreLOG2(hmm->ins[k].t[ts]);
 119           hmm->del[k].t[ts] = sreLOG2(hmm->del[k].t[ts]);
 120         }
 121       for (x = 0; x < Alphabet_size; x++)
 122         {
 123           hmm->mat[k].p[x]   = sreLOG2(hmm->mat[k].p[x]) - sreLOG2(hmm->null[x]);
 124           hmm->ins[k].p[x]   = sreLOG2(hmm->ins[k].p[x]) - sreLOG2(hmm->null[x]);
 125         }
 126     }
 127 }
 128
 129
 130 /* Function: WriteFlatPriorHMM()
 131  *
 132  * Purpose:  Fill an HMM with expected probabilities according
 133  *           to a given prior. Used to construct "flat" initial
 134  *           models for hmmt.
 135  */
 136 int
 137 WriteFlatPriorHMM(struct hmm_struc *hmm, struct prior_s *prior)
 138 {
 139   int    k;                     /* counter across model                */
 140   int    q;                     /* counter over mixtures               */
 141   int    x;                     /* counter over symbols or transitions */
 142   float  malpha;                /* alpha for mixture                   */
 143   float  ialpha;                /* alpha for insert mixture            */
 144   float  dalpha;                /* alpha for delete mixture            */
 145
 146   for (k = 0; k <= hmm->M; k++)
 147     {
 148                                 /* xray info for structure prior */
 149       if (prior->strategy == PRI_STRUCT)
 150         {
 151           hmm->xray[k*NINPUTS + XRAY_bias] = 1.0;
 152           hmm->xray[k*NINPUTS + XRAY_E]    = 0.0;
 153           hmm->xray[k*NINPUTS + XRAY_H]    = 0.0;
 154           hmm->xray[k*NINPUTS + XRAY_SA]   = 0.0;
 155         }
 156                                 /* match symbol emissions */
 157       for (x = 0; x < Alphabet_size; x++)
 158         hmm->mat[k].p[x] = 0.0;
 159       if (k > 0)
 160         for (q = 0; q < prior->mnum; q++)
 161           {
 162             if (prior->strategy == PRI_STRUCT)
 163               prior->mq[q] = 1.0 / prior->mnum;
 164             malpha = 0.0;
 165             for (x = 0; x < Alphabet_size; x++)
 166               malpha += prior->mat[q][x];
 167             for (x = 0; x < Alphabet_size; x++)
 168               hmm->mat[k].p[x] += prior->mq[q] * prior->mat[q][x] / malpha;
 169           }
 170                                 /* insert emissions */
 171       for (x = 0; x < Alphabet_size; x++)
 172         hmm->ins[k].p[x] = 0.0;
 173       for (q = 0; q < prior->inum; q++)
 174         {
 175           if (prior->strategy == PRI_STRUCT)
 176             prior->iq[q] = 1.0 / prior->inum;
 177           ialpha = 0.0;
 178           for (x = 0; x < Alphabet_size; x++)
 179             ialpha += prior->ins[q][x];
 180           for (x = 0; x < Alphabet_size; x++)
 181             hmm->ins[k].p[x] += prior->iq[q] * prior->ins[q][x] / ialpha;
 182         }
 183
 184                                 /* state transitions  */
 185       for (x = 0; x < 3; x++)
 186         hmm->mat[k].t[x] = hmm->ins[k].t[x] = hmm->del[k].t[x] = 0.0;
 187       for (q = 0; q < prior->tnum; q++)
 188         {
 189           if (prior->strategy == PRI_STRUCT)
 190             prior->tq[q] = 1.0 / prior->tnum;
 191           malpha = ialpha = dalpha = 0.0;
 192           for (x = 0; x < 3; x++)
 193             {
 194               malpha += prior->tm[q][x];
 195               ialpha += prior->ti[q][x];
 196               dalpha += prior->td[q][x];
 197             }
 198           for (x = 0; x < 3; x++)
 199             {
 200               hmm->mat[k].t[x] += prior->tq[q] * prior->tm[q][x] / malpha;
 201               hmm->ins[k].t[x] += prior->tq[q] * prior->ti[q][x] / ialpha;
 202               if (k > 0) hmm->del[k].t[x] += prior->tq[q] * prior->td[q][x] / dalpha;
 203             }
 204         }
 205     }
 206                                 /* the final state never transits to d+1 */
 207   hmm->mat[hmm->M].t[DELETE] = 0.0;
 208   hmm->ins[hmm->M].t[DELETE] = 0.0;
 209   hmm->del[hmm->M].t[DELETE] = 0.0;
 210   Renormalize(hmm);
 211   return 1;
 212 }
 213
 214
 215 /* Function: HMMDup()
 216  *
 217  * Purpose:  Create a duplicate copy of an HMM.
 218  *
 219  * Return:   Pointer to the duplicate.
 220  *           Caller is responsible for free'ing the duplicate.
 221  */
 222 struct hmm_struc *
 223 HMMDup(struct hmm_struc *hmm)
 224 {
 225   struct hmm_struc *newhmm;
 226
 227   if ((newhmm = AllocHMM(hmm->M)) == NULL)
 228     Die("AllocHMM() failed");
 229   HMMCopy(newhmm, hmm);
 230   return newhmm;
 231 }
 232
 233
 234 /* Function: HMMCopy()
 235  *
 236  * Purpose:  Make a copy of hmm2 in hmm1.
 237  *
 238  * Return:   (void)
 239  *           Caller promises that hmm1 and hmm2 have identical architectures.
 240  */
 241 void
 242 HMMCopy(struct hmm_struc *hmm1, struct hmm_struc *hmm2)
 243 {
 244   int k, x, ts;
 245
 246   hmm1->flags = hmm2->flags;
 247   if (hmm1->name != NULL) free(hmm1->name);
 248   hmm1->name = Strdup(hmm2->name);
 249
 250   if (hmm2->flags & HMM_REF)  strcpy(hmm1->ref, hmm2->ref);
 251   if (hmm2->flags & HMM_CS)   strcpy(hmm1->cs,  hmm2->cs);
 252   if (hmm2->flags & HMM_XRAY)
 253     memcpy(hmm1->xray, hmm2->xray, NINPUTS * (hmm2->M+2) * sizeof(float));
 254   memcpy(hmm1->null, hmm2->null, sizeof(float) * Alphabet_size);
 255
 256   for (k = 0; k <= hmm2->M+1; k++)
 257     {
 258                         /* copy transition T's */
 259       for (ts = 0; ts < 3; ts++)
 260         {
 261           hmm1->mat[k].t[ts] = hmm2->mat[k].t[ts];
 262           hmm1->ins[k].t[ts] = hmm2->ins[k].t[ts];
 263           hmm1->del[k].t[ts] = hmm2->del[k].t[ts];
 264         }
 265                                 /* copy symbol P tables */
 266       for (x = 0; x < Alphabet_size; x++)
 267         {
 268           hmm1->mat[k].p[x]   = hmm2->mat[k].p[x];
 269           hmm1->ins[k].p[x]   = hmm2->ins[k].p[x];
 270         }
 271     }
 272   return;
 273 }
 274
 275
 276 int
 277 FreeHMM(struct hmm_struc *hmm)
 278 {
 279   if (hmm == NULL) return 0;
 280   free(hmm->ref);
 281   free(hmm->cs);
 282   free(hmm->xray);
 283   free(hmm->name);
 284   if (hmm->mat != NULL)  free (hmm->mat);
 285   if (hmm->ins != NULL)  free (hmm->ins);
 286   if (hmm->del != NULL)  free (hmm->del);
 287   free(hmm);
 288   return 1;
 289 }
 290
 291
 292 struct shmm_s *
 293 AllocSearchHMM(int M)
 294 {
 295   struct shmm_s *shmm;
 296   int            x;
 297
 298   if ((shmm = (struct shmm_s *) malloc (sizeof(struct shmm_s))) == NULL)
 299     Die("malloc failed");
 300   for (x = 0; x < 26; x++)
 301     if ((shmm->m_emit[x] = (int *) calloc (M+1, sizeof(int))) == NULL ||
 302         (shmm->i_emit[x] = (int *) calloc (M+1, sizeof(int))) == NULL)
 303       Die("malloc failed");
 304   if ((shmm->t   = (int *)  malloc (sizeof(int)  * (9*(M+1)))) == NULL ||
 305       (shmm->ref = (char *) malloc (sizeof(char) * (M+2))) == NULL ||
 306       (shmm->cs  = (char *) malloc (sizeof(char) * (M+2))) == NULL)
 307     Die("malloc failed");
 308   shmm->flags = 0;
 309   shmm->name = Strdup("nameless");
 310   shmm->M = M;
 311   return shmm;
 312 }
 313
 314 void
 315 FreeSearchHMM(struct shmm_s *shmm)
 316 {
 317   int x;
 318
 319   for (x = 0; x < 26; x++)
 320     {
 321       free(shmm->m_emit[x]);
 322       free(shmm->i_emit[x]);
 323     }
 324   free(shmm->t);
 325   free(shmm->ref);
 326   free(shmm->cs);
 327   free(shmm->name);
 328   free(shmm);
 329 }
 330
 331
 332 /* Function: CountSymbol()
 333  *
 334  * Purpose:  Given an observed symbol, and a number of counts to
 335  *           distribute (typically just 1.0), bump the appropriate counter(s).
 336  *
 337  *           This is completely trivial only so long as the symbols
 338  *           always come from the expected alphabet; since we also
 339  *           have to deal with degenerate symbols for both nucleic
 340  *           acid and protein languages, we make a function to deal
 341  *           with this.
 342  *
 343  * Args:     sym      - observed symbol, e.g. `A' or `X'
 344  *           wt       - number of counts to distribute (e.g. 1.0)
 345  *           counters - array of 4 or 20 counters to increment
 346  *
 347  * Return:   Returns 1 on success and bumps the necessary counters.
 348  *           Returns 0 on failure and bumps each counter evenly, as
 349  *           if it saw a completely ambiguous symbol; this lets
 350  *           the caller silently accept garbage symbols, if it cares to.
 351  */
 352 int
 353 CountSymbol(char sym, float wt, float *counters)
 354 {
 355   char *sptr;                   /* pointer into symbol in hmm->alphabet         */
 356   int   status;                 /* RETURN: status; did we recognize the symbol? */
 357   char  symidx;                 /* index of symbol in Alphabet_iupac */
 358
 359   if ((sptr = strchr(Alphabet,sym)) != NULL)
 360     {
 361       symidx = (char) (sptr - Alphabet);
 362       status = 1;
 363     }
 364   else
 365     {
 366       symidx = (char) (Alphabet_iupac - 1);
 367       Warn("unrecognized character %c in CountSymbol()\n", sym);
 368       status = 0;
 369     }
 370   P7CountSymbol(counters, symidx, wt);
 371   return status;
 372 }
 373
 374
 375 /* Function: HMMDistance()
 376  *
 377  * Purpose:  Test two models for how different they are, using
 378  *           a simple squared difference measure on all homologous
 379  *           parameters. They must have the same architecture:
 380  *           i.e. check that newhmm->M == oldhmm->M before calling.
 381  *
 382  * Args:     newhmm  - new HMM, probability form
 383  *           oldhmm  - old HMM, probability form
 384  *
 385  * Return:   distance.
 386  */
 387 float
 388 HMMDistance(struct hmm_struc *newhmm, struct hmm_struc *oldhmm)
 389 {
 390   int    k,x, ts;
 391   float  distance = 0.0;
 392
 393   for (k = 0; k <= newhmm->M; k++)
 394     {
 395                                 /* state transition distances */
 396       if (k > 0)
 397         {
 398           for (ts = 0; ts < 3; ts++)
 399             distance += SQR( 100. * (newhmm->del[k].t[ts] - oldhmm->del[k].t[ts]));
 400         }
 401       for (ts = 0; ts < 3; ts++)
 402         distance += SQR( 100. * (newhmm->mat[k].t[ts] - oldhmm->mat[k].t[ts]));
 403       for (ts = 0; ts < 3; ts++)
 404         distance += SQR( 100. * (newhmm->ins[k].t[ts] - oldhmm->ins[k].t[ts]));
 405
 406                                 /* symbol emission distances */
 407       if (k > 0)
 408         for (x = 0; x < Alphabet_size; x++)
 409           distance += SQR( 100. * (newhmm->mat[k].p[x] - oldhmm->mat[k].p[x]));
 410       for (x = 0; x < Alphabet_size; x++)
 411           distance += SQR( 100. * (newhmm->ins[k].p[x] - oldhmm->ins[k].p[x]));
 412     }
 413   distance = sqrt(distance) / newhmm->M;
 414   return distance;
 415 }
 416
 417
 418
 419
 420 /* Function: Renormalize()
 421  *
 422  * Normalize all P distributions so they sum to 1.
 423  * P distributions that are all 0, or contain negative
 424  * probabilities, are left untouched.
 425  *
 426  * Returns 1 on success, or 0 on failure.
 427  */
 428 void
 429 Renormalize(struct hmm_struc *hmm)
 430 {
 431   int    k;                     /* counter for states                  */
 432
 433   for (k = 0; k <= hmm->M ; k++)
 434     {
 435                                 /* match state transition frequencies */
 436       FNorm(hmm->mat[k].t, 3);
 437       FNorm(hmm->ins[k].t, 3);
 438       if (k > 0) FNorm(hmm->del[k].t, 3);
 439
 440       if (k > 0) FNorm(hmm->mat[k].p, Alphabet_size);
 441       FNorm(hmm->ins[k].p, Alphabet_size);
 442     }
 443 }
 444