forester/archive/RIO/others/hmmer/src/debug.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 /* debug.c
  12  * Thu Nov 21 09:58:05 1996
  13  *
  14  * Printing out or naming various useful things from HMMER
  15  * innards.
  16  *
  17  * CVS $Id: debug.c,v 1.1.1.1 2005/03/22 08:33:58 cmzmasek Exp $
  18  */
  19
  20 #include <stdio.h>
  21 #include <stdlib.h>
  22 #include <stdarg.h>
  23 #include <ctype.h>
  24
  25 #include "structs.h"
  26 #include "config.h"
  27 #include "funcs.h"
  28 #include "squid.h"
  29
  30 /* Function: Statetype()
  31  *
  32  * Purpose:  Returns the state type in text.
  33  * Example:  Statetype(S) = "S"
  34  */
  35 char *
  36 Statetype(char st)
  37 {
  38   switch (st) {
  39   case STS: return "S";
  40   case STN: return "N";
  41   case STB: return "B";
  42   case STM: return "M";
  43   case STD: return "D";
  44   case STI: return "I";
  45   case STE: return "E";
  46   case STJ: return "J";
  47   case STC: return "C";
  48   case STT: return "T";
  49   default: return "BOGUS";
  50   }
  51 }
  52
  53 /* Function: AlphabetType2String()
  54  * Date:     SRE, Sun Dec 24 11:33:40 2000 [St. Louis]
  55  *
  56  * Purpose:  Returns a string "protein" for hmmAMINO,
  57  *           "nucleic acid" for hmmNUCLEIC, etc... used
  58  *           for formatting diagnostics.
  59  *
  60  * Args:     type - Alphabet type, e.g. hmmAMINO
  61  *
  62  * Returns:  char *
  63  */
  64 char *
  65 AlphabetType2String(int type)
  66 {
  67   switch (type) {
  68   case hmmAMINO:     return "protein";
  69   case hmmNUCLEIC:   return "nucleic acid";
  70   case hmmNOTSETYET: return "unknown";
  71   default:           return "BOGUS";
  72   }
  73 }
  74
  75
  76 /* Function: P7PrintTrace()
  77  *
  78  * Purpose:  Print out a traceback structure.
  79  *           If hmm is non-NULL, also print transition and emission scores.
  80  *
  81  * Args:     fp  - stderr or stdout, often
  82  *           tr  - trace structure to print
  83  *           hmm - NULL or hmm containing scores to print
  84  *           dsq - NULL or digitized sequence trace refers to.
  85  */
  86 void
  87 P7PrintTrace(FILE *fp, struct p7trace_s *tr, struct plan7_s *hmm, char *dsq)
  88 {
  89   int tpos;                     /* counter for trace position */
  90   int sym;
  91   int sc;
  92
  93   if (hmm == NULL) {
  94     fprintf(fp, "st  node   rpos  - traceback len %d\n", tr->tlen);
  95     fprintf(fp, "--  ---- ------\n");
  96     for (tpos = 0; tpos < tr->tlen; tpos++) {
  97       fprintf(fp, "%1s  %4d %6d\n",
  98               Statetype(tr->statetype[tpos]),
  99               tr->nodeidx[tpos],
 100               tr->pos[tpos]);
 101     }
 102   } else {
 103     if (!(hmm->flags & PLAN7_HASBITS))
 104       Die("oi, you can't print scores from that hmm, it's not ready.");
 105
 106     sc = 0;
 107     fprintf(fp, "st  node   rpos  transit emission - traceback len %d\n", tr->tlen);
 108     fprintf(fp, "--  ---- ------  ------- --------\n");
 109     for (tpos = 0; tpos < tr->tlen; tpos++) {
 110       if (dsq != NULL) sym = (int) dsq[tr->pos[tpos]];
 111
 112       fprintf(fp, "%1s  %4d %6d  %7d",
 113               Statetype(tr->statetype[tpos]),
 114               tr->nodeidx[tpos],
 115               tr->pos[tpos],
 116               (tpos < tr->tlen-1) ?
 117               TransitionScoreLookup(hmm, tr->statetype[tpos], tr->nodeidx[tpos],
 118                                     tr->statetype[tpos+1], tr->nodeidx[tpos+1]) : 0);
 119
 120       if (tpos < tr->tlen-1)
 121         sc += TransitionScoreLookup(hmm, tr->statetype[tpos], tr->nodeidx[tpos],
 122                                     tr->statetype[tpos+1], tr->nodeidx[tpos+1]);
 123
 124       if (dsq != NULL) {
 125         if (tr->statetype[tpos] == STM)
 126           {
 127             fprintf(fp, " %8d %c", hmm->msc[sym][tr->nodeidx[tpos]],
 128                     Alphabet[sym]);
 129             sc += hmm->msc[sym][tr->nodeidx[tpos]];
 130           }
 131         else if (tr->statetype[tpos] == STI)
 132           {
 133             fprintf(fp, " %8d %c", hmm->isc[sym][tr->nodeidx[tpos]],
 134                     (char) tolower((int) Alphabet[sym]));
 135             sc += hmm->isc[sym][tr->nodeidx[tpos]];
 136           }
 137         else if ((tr->statetype[tpos] == STN && tr->statetype[tpos-1] == STN) ||
 138                  (tr->statetype[tpos] == STC && tr->statetype[tpos-1] == STC) ||
 139                  (tr->statetype[tpos] == STJ && tr->statetype[tpos-1] == STJ))
 140           {
 141             fprintf(fp, " %8d %c", 0, (char) tolower((int) Alphabet[sym]));
 142           }
 143       } else {
 144         fprintf(fp, " %8s %c", "-", '-');
 145       }
 146
 147
 148       fputs("\n", fp);
 149     }
 150     fprintf(fp, "                 ------- --------\n");
 151     fprintf(fp, "           total: %6d\n\n", sc);
 152   }
 153 }
 154
 155 /* Function: P7PrintPrior()
 156  *
 157  * Purpose:  Print out a Plan 7 prior structure.
 158  */
 159 void
 160 P7PrintPrior(FILE *fp, struct p7prior_s *pri)
 161 {
 162   int q, x;                     /* counters for mixture component, element */
 163
 164   if      (pri->strategy == PRI_DCHLET) fputs("Dirichlet\n", fp);
 165   else if (pri->strategy == PRI_PAM)    fputs("PAM\n", fp);
 166   else Die("No such strategy.");
 167
 168   if      (Alphabet_type == hmmAMINO)   fputs("Amino\n", fp);
 169   else if (Alphabet_type == hmmNUCLEIC) fputs("Nucleic\n", fp);
 170
 171   /* Transitions
 172    */
 173   fprintf(fp, "\n%d\n", pri->tnum);
 174   for (q = 0; q < pri->tnum; q++)
 175     {
 176       fprintf(fp, "%.4f\n", pri->tq[q]);
 177       for (x = 0; x < 7; x++)
 178         fprintf(fp, "%.4f ", pri->t[q][x]);
 179       fputs("\n", fp);
 180     }
 181
 182   /* Match emissions
 183    */
 184   fprintf(fp, "\n%d\n", pri->mnum);
 185   for (q = 0; q < pri->mnum; q++)
 186     {
 187       fprintf(fp, "%.4f\n", pri->mq[q]);
 188       for (x = 0; x < Alphabet_size; x++)
 189         fprintf(fp, "%.4f ", pri->m[q][x]);
 190       fputs("\n", fp);
 191     }
 192
 193   /* Insert emissions
 194    */
 195   fprintf(fp, "\n%d\n", pri->inum);
 196   for (q = 0; q < pri->inum; q++)
 197     {
 198       fprintf(fp, "%.4f\n", pri->iq[q]);
 199       for (x = 0; x < Alphabet_size; x++)
 200         fprintf(fp, "%.4f ", pri->i[q][x]);
 201       fputs("\n", fp);
 202     }
 203 }
 204
 205 /* Function: TraceVerify()
 206  * Date:     SRE, Mon Feb  2 07:48:52 1998 [St. Louis]
 207  *
 208  * Purpose:  Check a traceback structure for internal consistency.
 209  *           Used in Shiva testsuite, for example.
 210  *
 211  * Args:     tr  - traceback to verify
 212  *           M   - length of HMM
 213  *           N   - length of sequence
 214  *
 215  * Returns:  1 if OK. 0 if not.
 216  */
 217 int
 218 TraceVerify(struct p7trace_s *tr, int M, int N)
 219 {
 220   int tpos;                     /* position in trace                  */
 221   int k;                        /* current position in HMM nodes 1..M */
 222   int i;                        /* current position in seq 1..N       */
 223   int nn, nc, nj;               /* number of STN's, STC's, STJ's seen */
 224   int nm;                       /* number of STM's seen */
 225
 226   /* Basic checks on ends.
 227    */
 228   if (tr->statetype[0] != STS)          return 0;
 229   if (tr->statetype[1] != STN)          return 0;
 230   if (tr->statetype[tr->tlen-2] != STC) return 0;
 231   if (tr->statetype[tr->tlen-1] != STT) return 0;
 232   if (tr->pos[1] != 0)                  return 0;
 233
 234   /* Check for consistency throughout trace
 235    */
 236   k = i = nn = nc = nj = nm = 0;
 237   for (tpos = 0; tpos < tr->tlen; tpos++)
 238     {
 239       switch (tr->statetype[tpos]) {
 240       case STS:
 241         if (tr->nodeidx[tpos] != 0) return 0;
 242         if (tr->pos[tpos]     != 0) return 0;
 243         if (k != 0)                 return 0;
 244         if (i != 0)                 return 0;
 245         if (tpos != 0)              return 0;
 246         break;
 247
 248       case STN:                 /* first N doesn't emit. */
 249         if (tr->nodeidx[tpos] != 0) return 0;
 250         if (k != 0)                 return 0;
 251         if (nn > 0)
 252           {
 253             if (tr->pos[tpos] != i+1) return 0;
 254             i++;
 255           }
 256         else
 257           {
 258             if (tr->pos[tpos] != 0) return 0;
 259             if (i != 0)             return 0;
 260           }
 261         nn++;
 262         break;
 263
 264       case STB:
 265         if (tr->nodeidx[tpos] != 0) return 0;
 266         if (tr->pos[tpos]     != 0) return 0;
 267         nm = 0;
 268         break;
 269
 270       case STM:                 /* can enter anywhere on first M */
 271         if (tr->pos[tpos] != i+1) return 0;
 272         if (tr->nodeidx[tpos] < 1 || tr->nodeidx[tpos] > M) return 0;
 273         i++;
 274         if (nm == 0)  k = tr->nodeidx[tpos];
 275         else {
 276           if (tr->nodeidx[tpos] != k+1) return 0;
 277           k++;
 278         }
 279         nm++;
 280         break;
 281
 282       case STI:
 283         if (tr->pos[tpos] != i+1)   return 0;
 284         if (tr->nodeidx[tpos] != k) return 0;
 285         if (tr->nodeidx[tpos] < 1 || tr->nodeidx[tpos] > M-1) return 0;
 286         if (k >= M)                 return 0;
 287         i++;
 288         break;
 289
 290       case STD:
 291         if (tr->pos[tpos] != 0)       return 0;
 292         if (tr->nodeidx[tpos] != k+1) return 0;
 293         if (tr->nodeidx[tpos] < 1 || tr->nodeidx[tpos] > M) return 0;
 294         k++;
 295         break;
 296
 297       case STE:
 298         if (tr->nodeidx[tpos] != 0) return 0;
 299         if (tr->pos[tpos]     != 0) return 0;
 300         nj = 0;
 301         break;
 302
 303       case STJ:
 304         if (tr->nodeidx[tpos] != 0) return 0;
 305         if (nj > 0)
 306           {
 307             if (tr->pos[tpos] != i+1) return 0;
 308             i++;
 309           }
 310         else if (tr->pos[tpos] != 0) return 0;
 311         nj++;
 312         break;
 313
 314       case STC:
 315         if (tr->nodeidx[tpos] != 0) return 0;
 316         if (nc > 0)
 317           {
 318             if (tr->pos[tpos] != i+1) return 0;
 319             i++;
 320           }
 321         else if (tr->pos[tpos] != 0)  return 0;
 322         nc++;
 323         break;
 324
 325       case STT:
 326         if (tpos != tr->tlen - 1)   return 0;
 327         if (tr->nodeidx[tpos] != 0) return 0;
 328         if (tr->pos[tpos]     != 0) return 0;
 329         if (i != N)                 return 0;
 330         break;
 331
 332       case STBOGUS:
 333       default:
 334         return 0;
 335       } /* end switch over statetypes */
 336     } /* end loop over trace positions */
 337
 338   return 1;
 339 }
 340
 341
 342 /* Function: TraceCompare()
 343  * Date:     SRE, Wed Mar  4 17:26:49 1998 [St. Louis]
 344  *
 345  * Purpose:  Compare two tracebacks; return 1 if they're
 346  *           identical, else 0. Written for Shiva testsuite.
 347  *
 348  * Args:     t1 - first trace
 349  *           t2 - second trace
 350  *
 351  * Returns:  1 if identical; 0 elsewise
 352  */
 353 int
 354 TraceCompare(struct p7trace_s *t1, struct p7trace_s *t2)
 355 {
 356   int tpos;
 357
 358   if (t1->tlen != t2->tlen) return 0;
 359
 360   for (tpos = 0; tpos < t1->tlen; tpos++)
 361     {
 362       if (t1->statetype[tpos] != t2->statetype[tpos]) return 0;
 363       if (t1->nodeidx[tpos]   != t2->nodeidx[tpos])   return 0;
 364       if (t1->pos[tpos]       != t2->pos[tpos])       return 0;
 365     }
 366   return 1;
 367 }
 368