forester/archive/RIO/others/hmmer/testsuite/trace_test.c

   1 /* trace_test.c
   2  * Mon Feb  2 07:57:47 1998
   3  * cp trace_test.c ../src/testdriver.c; cd ../src; make testdriver
   4  *
   5  * Test driver for Viterbi tracebacks.
   6  *
   7  * RCS $Id: trace_test.c,v 1.1.1.1 2005/03/22 08:34:47 cmzmasek Exp $
   8  */
   9
  10
  11 #include <stdio.h>
  12 #include <time.h>
  13 #include <math.h>
  14
  15 #include "structs.h"
  16 #include "funcs.h"
  17 #include "globals.h"
  18 #include "squid.h"
  19
  20 static char banner[] = "\
  21 trace_test : testing of Plan7 Viterbi traceback code";
  22
  23 static char usage[] = "\
  24 Usage: testdriver [-options]\n\
  25   Available options are:\n\
  26   -h              : help; display this usage info\n\
  27   -v              : be verbose\n\
  28 ";
  29
  30 static char experts[] = "\
  31   --hmm <f>       : use HMM in file <f>\n\
  32   --seq <f>       : use seq(s) in file <f>\n\
  33   --small         : run P7SmallViterbi()\n\
  34 \n";
  35
  36 static struct opt_s OPTIONS[] = {
  37   { "-h",       TRUE,  sqdARG_NONE },
  38   { "-v",       TRUE,  sqdARG_NONE },
  39   { "--hmm",    FALSE, sqdARG_STRING },
  40   { "--seq",    FALSE, sqdARG_STRING },
  41   { "--small",  FALSE, sqdARG_NONE },
  42 };
  43 #define NOPTIONS (sizeof(OPTIONS) / sizeof(struct opt_s))
  44
  45 int
  46 main(int argc, char **argv)
  47 {
  48   char    *hmmfile;             /* file to read HMM(s) from                */
  49   HMMFILE *hmmfp;               /* opened hmmfile for reading              */
  50   char    *seqfile;             /* file to read target sequence(s) from    */
  51   SQFILE   *sqfp;               /* opened seqfile for reading              */
  52   char     *seq;                /* target sequence                         */
  53   SQINFO    sqinfo;             /* optional info for seq                   */
  54   char     *dsq;                /* digitized target sequence               */
  55   struct plan7_s  *hmm;         /* HMM to search with                      */
  56   struct p7trace_s  *tr;        /* traceback                               */
  57   int       nseq;
  58   float     sc;
  59
  60   int be_verbose;
  61   int do_small;                 /* TRUE to invoke P7SmallViterbi */
  62
  63   char *optname;                /* name of option found by Getopt()         */
  64   char *optarg;                 /* argument found by Getopt()               */
  65   int   optind;                 /* index in argv[]                          */
  66
  67   /***********************************************
  68    * Parse command line
  69    ***********************************************/
  70
  71   be_verbose = FALSE;
  72   hmmfile    = "trace_test.hmm";
  73   seqfile    = "trace_test.seq";
  74   do_small   = FALSE;
  75
  76   while (Getopt(argc, argv, OPTIONS, NOPTIONS, usage,
  77                 &optind, &optname, &optarg))  {
  78     if      (strcmp(optname, "-v")       == 0) be_verbose = TRUE;
  79     else if (strcmp(optname, "--hmm")    == 0) hmmfile    = optarg;
  80     else if (strcmp(optname, "--seq")    == 0) seqfile    = optarg;
  81     else if (strcmp(optname, "--small")  == 0) do_small   = TRUE;
  82     else if (strcmp(optname, "-h")       == 0) {
  83       Banner(stdout, banner);
  84       puts(usage);
  85       puts(experts);
  86       exit(0);
  87     }
  88   }
  89   if (argc - optind != 0)
  90     Die("Incorrect number of arguments.\n%s\n", usage);
  91
  92   /***********************************************
  93    * Open test sequence file
  94    ***********************************************/
  95
  96   if ((sqfp = SeqfileOpen(seqfile, SQFILE_UNKNOWN, "BLASTDB")) == NULL)
  97     Die("Failed to open sequence database file %s\n%s\n", seqfile, usage);
  98
  99   /***********************************************
 100    * Open HMM file
 101    * Read a single HMM from it. (Config HMM, if necessary).
 102    ***********************************************/
 103
 104   if ((hmmfp = HMMFileOpen(hmmfile, NULL)) == NULL)
 105     Die("Failed to open HMM file %s\n%s", hmmfile, usage);
 106   if (!HMMFileRead(hmmfp, &hmm))
 107     Die("Failed to read any HMMs from %s\n", hmmfile);
 108   if (hmm == NULL)
 109     Die("HMM file %s corrupt or in incorrect format? Parse failed", hmmfile);
 110   P7Logoddsify(hmm, TRUE);
 111
 112   /***********************************************
 113    * Search HMM against each sequence
 114    ***********************************************/
 115
 116   nseq = 0;
 117   while (ReadSeq(sqfp, sqfp->format, &seq, &sqinfo))
 118     {
 119       nseq++;
 120       dsq = DigitizeSequence(seq, sqinfo.len);
 121
 122       if (do_small) sc = P7SmallViterbi(dsq, sqinfo.len, hmm, &tr);
 123       else          sc = P7Viterbi(dsq, sqinfo.len, hmm, &tr);
 124
 125       if (be_verbose)
 126         {
 127           printf("test sequence %d: score %.1f : %s %s\n",
 128                  nseq, sc, sqinfo.name,
 129                  sqinfo.flags & SQINFO_DESC ? sqinfo.desc : "");
 130           P7PrintTrace(stdout, tr, hmm, dsq);
 131         }
 132
 133       if (! TraceVerify(tr, hmm->M, sqinfo.len))
 134         Die("Trace verify failed on seq #%d, %s\n", nseq, sqinfo.name);
 135
 136       FreeSequence(seq, &sqinfo);
 137       P7FreeTrace(tr);
 138       free(dsq);
 139     }
 140
 141   FreePlan7(hmm);
 142   HMMFileClose(hmmfp);
 143   SeqfileClose(sqfp);
 144
 145   return EXIT_SUCCESS;
 146 }