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

   1 /************************************************************
   2  * Copyright (C) 1998 Ian Holmes
   3  * HMMER - Biological sequence analysis with profile HMMs
   4  * Copyright (C) 1992-1999 Washington University School of Medicine
   5  * All Rights Reserved
   6  *
   7  *     This source code is distributed under the terms of the
   8  *     GNU General Public License. See the files COPYING and LICENSE
   9  *     for details.
  10  ************************************************************/
  11
  12 /* display.c
  13  * Author: Ian Holmes (ihh@sanger.ac.uk, Jun 5 1998)
  14  * Derived from core_algorithms.c (SRE, Nov 11 1996)
  15  * Incorporated SRE, Sat Nov  6 10:09:41 1999
  16  *
  17  * Functions for displaying HMMer2.0 structures.
  18  *
  19  * RCS $Id: display.c,v 1.1.1.1 2005/03/22 08:33:59 cmzmasek Exp $
  20  */
  21
  22 #include "structs.h"
  23 #include "config.h"
  24 #include "funcs.h"
  25 #include "squid.h"
  26
  27 void PrintIscore(int sc);
  28
  29 void PrintTransition(char src,
  30                      int isrc,
  31                      int ksrc,
  32                      char dest,
  33                      int idest,
  34                      int kdest,
  35                      int sc,
  36                      struct p7trace_s **alignment,
  37                      int *min,
  38                      int *max,
  39                      int *on,
  40                      int A);
  41
  42
  43 /* Function: DisplayPlan7Posteriors()
  44  *
  45  * Purpose:  Print out posterior transition probabilities
  46  *           in modelpost format.
  47  *           NB only prints out transitions that touch
  48  *           either the Viterbi or the optimal accuracy path.
  49  *
  50  * Args:     L        - the length of the sequence
  51  *           hmm      - the model
  52  *           forward  - forward matrix
  53  *           backward - backward matrix
  54  *           viterbi  - Viterbi trace
  55  *           optacc   - optimal accuracy trace
  56  *
  57  * Return:   void
  58  *
  59  */
  60 void DisplayPlan7Posteriors(int L, struct plan7_s *hmm,
  61                             struct dpmatrix_s *forward,
  62                             struct dpmatrix_s *backward,
  63                             struct p7trace_s *viterbi,
  64                             struct p7trace_s *optacc)
  65 {
  66   struct p7trace_s* alignment[2];
  67   alignment[0] = viterbi;
  68   alignment[1] = optacc;
  69   DisplayPlan7PostAlign (L, hmm, forward, backward, alignment, 2);
  70 }
  71
  72
  73 /* Function: DisplayPlan7PostAlign()
  74  *
  75  * Purpose:  Print out posterior transition probabilities
  76  *           in modelpost format, for any set of alignments.
  77  *
  78  * Args:     L         - the length of the sequence
  79  *           hmm       - the model
  80  *           forward   - forward matrix
  81  *           backward  - backward matrix
  82  *           alignment - array of traces
  83  *           A         - size of alignment array
  84  *
  85  * Return:   void
  86  *
  87  */
  88 void DisplayPlan7PostAlign(int L, struct plan7_s *hmm,
  89                            struct dpmatrix_s *forward,
  90                            struct dpmatrix_s *backward,
  91                            struct p7trace_s **alignment,
  92                            int A)
  93 {
  94   int sc;
  95   int i;
  96   int j;
  97   int k;
  98   int kmin;
  99   int kmax;
 100   int* min;
 101   int* max;
 102   int* on;
 103   char state;
 104
 105   sc = forward->xmx[L][XMC] + hmm->xsc[XTC][MOVE];     /* total Forward score */
 106
 107   min = (int*) calloc (A, sizeof(int));
 108   max = (int*) calloc (A, sizeof(int));
 109   on  = (int*) calloc (A, sizeof(int));
 110
 111   for (i = 0; i <= L; i++)
 112     {
 113       for (j = 0; j < A; j++) {
 114         while (alignment[j]->pos[min[j]] < i - 1 && min[j] < alignment[j]->tlen - 1)
 115           min[j]++;
 116
 117         while (alignment[j]->pos[max[j]] <= i + 1 && max[j] < alignment[j]->tlen - 1)
 118           max[j]++;
 119       }
 120
 121       for (state = STM; state <= STJ; state++)
 122         {
 123           if (state == STM || state == STB)
 124             {
 125               kmin = 1;
 126               kmax = hmm->M;
 127             }
 128           else if (state == STD)
 129             {
 130               kmin = 2;
 131               kmax = hmm->M - 1;
 132             }
 133           else if (state == STI)
 134             {
 135               kmin = 1;
 136               kmax = hmm->M - 1;
 137             }
 138           else
 139             kmin = kmax = 0;
 140
 141           for (k = kmin; k <= kmax; k++)
 142             {
 143               switch (state)
 144                 {
 145                 case STM:
 146                   if (i<L && k<hmm->M)
 147                     PrintTransition (STM,i,k, STM,i+1,k+1,
 148                                      forward->mmx[i][k] + hmm->tsc[k][TMM] + backward->mmx[i+1][k+1] - sc,
 149                                      alignment, min, max, on, A);
 150
 151                   if (i<L && k<hmm->M)
 152                     PrintTransition (STM,i,k, STI,i+1,k,
 153                                      forward->mmx[i][k] + hmm->tsc[k][TMI] + backward->imx[i+1][k] - sc,
 154                                      alignment, min, max, on, A);
 155
 156                   if (k<hmm->M-1)
 157                     PrintTransition (STM,i,k, STD,i,k+1,
 158                                      forward->mmx[i][k] + hmm->tsc[k][TMD] + backward->dmx[i][k+1] - sc,
 159                                      alignment, min, max, on, A);
 160
 161                   PrintTransition (STM,i,k, STE,i,0,
 162                                    forward->mmx[i][k] + hmm->esc[k] + backward->xmx[i][XME] - sc,
 163                                    alignment, min, max, on, A);
 164                   break;
 165
 166                 case STD:
 167                   if (i<L)
 168                     PrintTransition (STD,i,k, STM,i+1,k+1,
 169                                      forward->dmx[i][k] + hmm->tsc[k][TDM] + backward->mmx[i+1][k+1] - sc,
 170                                      alignment, min, max, on, A);
 171
 172                   PrintTransition (STD,i,k, STD,i,k+1,
 173                                    forward->dmx[i][k] + hmm->tsc[k][TDD] + backward->dmx[i][k+1] - sc,
 174                                    alignment, min, max, on, A);
 175
 176                   break;
 177
 178                 case STI:
 179                   if (i<L)
 180                     PrintTransition (STI,i,k, STM,i+1,k+1,
 181                                      forward->imx[i][k] + hmm->tsc[k][TIM] + backward->mmx[i+1][k+1] - sc,
 182                                      alignment, min, max, on, A);
 183
 184                   if (i<L)
 185                     PrintTransition (STI,i,k, STI,i+1,k,
 186                                      forward->imx[i][k] + hmm->tsc[k][TII] + backward->imx[i+1][k] - sc,
 187                                      alignment, min, max, on, A);
 188
 189                   break;
 190
 191                 case STB:
 192                   if (i<L)
 193                     PrintTransition (STB,i,0, STM,i+1,k,
 194                                      forward->xmx[i][XMB] + hmm->bsc[k] + backward->mmx[i+1][k] - sc,
 195                                      alignment, min, max, on, A);
 196                   break;
 197
 198                 default:
 199                   break;
 200
 201                 }
 202             }
 203
 204           switch (state)
 205             {
 206             case STN:
 207               PrintTransition (STN,i,0, STB,i,0,
 208                                forward->xmx[i][XMN] + hmm->xsc[XTN][MOVE] + backward->xmx[i][XMB] - sc,
 209                                alignment, min, max, on, A);
 210
 211               if (i<L)
 212                 PrintTransition (STN,i,0, STN,i+1,0,
 213                                  forward->xmx[i][XMN] + hmm->xsc[XTN][LOOP] + backward->xmx[i+1][XMN] - sc,
 214                                  alignment, min, max, on, A);
 215               break;
 216
 217             case STJ:
 218               PrintTransition (STJ,i,0, STB,i,0,
 219                                forward->xmx[i][XMJ] + hmm->xsc[XTJ][MOVE] + backward->xmx[i][XMB] - sc,
 220                                alignment, min, max, on, A);
 221
 222               if (i<L)
 223                 PrintTransition (STJ,i,0, STJ,i+1,0,
 224                                  forward->xmx[i][XMJ] + hmm->xsc[XTJ][LOOP] + backward->xmx[i+1][XMJ] - sc,
 225                                  alignment, min, max, on, A);
 226               break;
 227
 228             case STC:
 229               PrintTransition (STC,i,0, STT,i,0,
 230                                forward->xmx[i][XMC] + hmm->xsc[XTC][MOVE] - sc,      /* should be 1 */
 231                                alignment, min, max, on, A);
 232
 233               if (i<L)
 234                 PrintTransition (STC,i,0, STC,i+1,0,
 235                                  forward->xmx[i][XMC] + hmm->xsc[XTC][LOOP] + backward->xmx[i+1][XMC] - sc,
 236                                  alignment, min, max, on, A);
 237               break;
 238
 239             case STE:
 240               PrintTransition (STE,i,0, STC,i,0,
 241                                forward->xmx[i][XME] + hmm->xsc[XTE][MOVE] + backward->xmx[i][XMC] - sc,
 242                                alignment, min, max, on, A);
 243
 244               PrintTransition (STE,i,0, STJ,i,0,
 245                                forward->xmx[i][XME] + hmm->xsc[XTE][LOOP] + backward->xmx[i][XMJ] - sc,
 246                                alignment, min, max, on, A);
 247               break;
 248
 249             case STS:
 250               if (i == 0)
 251                 PrintTransition (STS,i,0, STN,i,0,
 252                                  backward->xmx[i][XMN] - sc,          /* should be 1 */
 253                                  alignment, min, max, on, A);
 254               break;
 255
 256             case STM:
 257             case STD:
 258             case STI:
 259             case STB:
 260             case STT:
 261               break;
 262
 263             default:
 264               Die ("unknown state");
 265
 266             }
 267         }
 268     }
 269
 270   free (min);
 271   free (max);
 272   free (on);
 273
 274 }
 275
 276
 277
 278 /* Function: DisplayPlan7Matrix()
 279  *
 280  * Purpose:  Print out a dynamic programming matrix.
 281  *
 282  * Args:     dsq    - sequence in digitized form
 283  *           L      - length of dsq
 284  *           hmm    - the model
 285  *           mx     - dp matrix
 286  *
 287  * Return:   void
 288  *
 289  * The output of this function inverts HMMer's concept of rows and columns
 290  * (i.e. each row represents a state, and each column, a residue);
 291  * also, probabilities are displayed as natural logs, not bit scores.
 292  * It should probably only be used by ihh...
 293  *
 294  */
 295 void
 296 DisplayPlan7Matrix(char *dsq, int L, struct plan7_s *hmm, struct dpmatrix_s *mx)
 297 {
 298   int i;
 299   int k;
 300
 301   printf("         *      ");
 302   for (i=1;i<=L;i++) printf("    %c      ",Alphabet[dsq[i]]);
 303   printf("\nN    ");
 304   for (i=0;i<=L;i++) PrintIscore(mx->xmx[i][XMN]);
 305   for (k=1;k<=hmm->M;k++) {
 306     printf("\nM%-3d ",k);
 307     for (i=0;i<=L;i++) PrintIscore(mx->mmx[i][k]);
 308   }
 309   for (k=1;k<hmm->M;k++) {
 310     printf("\nI%-3d ",k);
 311     for (i=0;i<=L;i++) PrintIscore(mx->imx[i][k]);
 312   }
 313   printf("\nE    ");
 314   for (i=0;i<=L;i++) PrintIscore(mx->xmx[i][XME]);
 315   printf("\nC    ");
 316   for (i=0;i<=L;i++) PrintIscore(mx->xmx[i][XMC]);
 317   printf("\nJ    ");
 318   for (i=0;i<=L;i++) PrintIscore(mx->xmx[i][XMJ]);
 319   printf("\nB    ");
 320   for (i=0;i<=L;i++) PrintIscore(mx->xmx[i][XMB]);
 321   for (k=2;k<hmm->M;k++) {
 322     printf("\nD%-3d ",k);
 323     for (i=0;i<=L;i++) PrintIscore(mx->dmx[i][k]);
 324   }
 325   printf("\n\n");
 326 }
 327
 328
 329 void PrintIscore(int sc) {
 330   double dsc;
 331   double div;
 332   dsc = (double) sc;
 333   div = INTSCALE / 0.693147180559945;   /* == INTSCALE / log(2) */
 334   dsc = dsc / div;
 335   printf("%- #11.3e",dsc);
 336 }
 337
 338
 339 void PrintTransition(char src,
 340                      int isrc,
 341                      int ksrc,
 342                      char dest,
 343                      int idest,
 344                      int kdest,
 345                      int sc,
 346                      struct p7trace_s **alignment,
 347                      int *min,
 348                      int *max,
 349                      int *on,
 350                      int A)
 351 {
 352   char src_str[6];     /* buffer for source state label        */
 353   char dest_str[6];    /* buffer for destination state label   */
 354   int j;
 355   int tpos;
 356   int tnext;
 357   int pos;
 358   int next;
 359   int near;
 360
 361   near = 0;
 362
 363   for (j = 0; j < A; j++) {
 364     on[j] = 0;
 365     for (pos = 0, tpos = min[j]; tpos <= max[j]; tpos++) {
 366
 367       if (alignment[j]->pos[tpos] != 0)
 368         pos = alignment[j]->pos[tpos];
 369
 370       if (src == alignment[j]->statetype[tpos]
 371           && ksrc == alignment[j]->nodeidx[tpos]
 372           && isrc == pos)
 373         near = TRUE;
 374
 375       if (dest == alignment[j]->statetype[tpos]
 376           && kdest == alignment[j]->nodeidx[tpos]
 377           && idest == pos)
 378         near = TRUE;
 379
 380       if (tpos < alignment[j]->tlen - 1)
 381         {
 382           tnext = tpos + 1;
 383
 384           /* fold up B->D->M transitions into pseudo- B->M transitions */
 385
 386           if (alignment[j]->statetype[tpos] == STB)
 387             while (alignment[j]->statetype[tnext] == STD && tnext < alignment[j]->tlen - 1)
 388               tnext++;
 389
 390           next = alignment[j]->pos[tnext];
 391           if (next == 0)
 392             next = pos;
 393
 394           if (src == alignment[j]->statetype[tpos]
 395               && ksrc == alignment[j]->nodeidx[tpos]
 396               && isrc == pos
 397               && dest == alignment[j]->statetype[tnext]
 398               && kdest == alignment[j]->nodeidx[tnext]
 399               && idest == next)
 400             on[j] = TRUE;
 401         }
 402     }
 403   }
 404
 405   if (!near) return;
 406
 407   switch (src)
 408     {
 409     case STM: sprintf (src_str, "M%d", ksrc); break;
 410     case STD: sprintf (src_str, "D%d", ksrc); break;
 411     case STI: sprintf (src_str, "I%d", ksrc); break;
 412     case STS: sprintf (src_str, "S"); break;
 413     case STN: sprintf (src_str, "N"); break;
 414     case STB: sprintf (src_str, "B"); break;
 415     case STE: sprintf (src_str, "E"); break;
 416     case STC: sprintf (src_str, "C"); break;
 417     case STJ: sprintf (src_str, "J"); break;
 418     case STT: sprintf (src_str, "T"); break;
 419     default: Die ("bad transition");
 420     }
 421
 422   switch (dest)
 423     {
 424     case STM: sprintf (dest_str, "M%d", kdest); break;
 425     case STD: sprintf (dest_str, "D%d", kdest); break;
 426     case STI: sprintf (dest_str, "I%d", kdest); break;
 427     case STS: sprintf (dest_str, "S"); break;
 428     case STN: sprintf (dest_str, "N"); break;
 429     case STB: sprintf (dest_str, "B"); break;
 430     case STE: sprintf (dest_str, "E"); break;
 431     case STC: sprintf (dest_str, "C"); break;
 432     case STJ: sprintf (dest_str, "J"); break;
 433     case STT: sprintf (dest_str, "T"); break;
 434     default: Die ("bad transition");
 435     }
 436
 437   printf ("%d\t%s\t%d\t%s\t%-14.7g\t", isrc, src_str, idest, dest_str, (double) Score2Prob(sc,1.));
 438
 439   for (j = 0; j < A; j++) {
 440     if (on[j]) printf ("*");
 441     if (j < A - 1) printf ("\t");
 442   }
 443
 444   printf ("\n");
 445
 446 }
 447