forester/archive/RIO/others/hmmer/squid/rk.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 /* rk.c (originally from rnabob's patsearch.c)
  12  *
  13  * Contains a compiler and a search engine for Rabin-Karp
  14  * based primary sequence pattern searching on encoded
  15  * sequences.
  16  *
  17  * See Sedgewick, _Algorithms_, for a general discussion of
  18  * the Rabin-Karp algorithm. See the rkcomp or rkexec man
  19  * pages for specific details.
  20  *
  21  * RCS $Id: rk.c,v 1.1.1.1 2005/03/22 08:34:16 cmzmasek Exp $
  22  */
  23
  24 #include <stdio.h>
  25 #include <string.h>
  26 #include <ctype.h>
  27 #include "squid.h"              /* seq encoding utilities and typedefs */
  28 #include "rk.h"
  29
  30
  31 #ifdef MEMDEBUG
  32 #include "dbmalloc.h"
  33 #endif
  34
  35 Hashseq
  36 rkcomp(char *probe)               /* A,C,G,T/U, N probe string, 0-8 nt long */
  37 {
  38   Hashseq hashprobe = 0;
  39   char    coded[RK_HASHSIZE + 1];
  40   int     len;
  41   int     i;
  42                                 /* check bounds violation on probe */
  43   if ((len =  strlen(probe)) > RK_HASHSIZE) return 0;
  44                                 /* encode the probe */
  45   if (seqencode(coded, probe) == 0) return 0;
  46                                 /* pack the probe into a Hashseq */
  47   for (i = 0; i < len; i++)
  48     {
  49       hashprobe <<= 4;
  50       hashprobe |= (Hashseq) coded[i];
  51     }
  52                                 /* left adjust as needed */
  53   for (; i < RK_HASHSIZE; i++)
  54     {
  55       hashprobe <<= 4;
  56       hashprobe |= (Hashseq) NTN;
  57     }
  58                                 /* return the compiled probe */
  59   return hashprobe;
  60 }
  61
  62 int
  63 rkseq(Hashseq   hashprobe,      /* up to 8 nt packed into the probe */
  64       char     *sequence)       /* encoded sequence                 */
  65 {
  66   long     i;
  67   long     pos = 0;
  68   Hashseq  target = 0;
  69
  70                                 /* initialize the target hashseq */
  71   for (i = 0; i < RK_HASHSIZE; i++)
  72     {
  73       if (*(sequence + i) == NTEND)
  74         break;
  75       target <<= 4;
  76       target |=  (Hashseq) (*(sequence + i));
  77     }
  78
  79   while (*(sequence + pos + RK_HASHSIZE -1) != NTEND)
  80     {
  81 #ifdef DEBUG
  82       printf("hashprobe: ");
  83       writehash(hashprobe);
  84       printf("\ttarget: ");
  85       writehash(target);
  86       printf("\nhashprobe & target: ");
  87       writehash(hashprobe & target);
  88       printf("\n");
  89 #endif
  90       if ((hashprobe & target) == target)
  91         return ((int) pos);
  92       target <<= 4;
  93       target |=  (Hashseq) (*(sequence + pos + RK_HASHSIZE));
  94       pos++;
  95     }
  96                                 /* now we deal with an end effect */
  97   for (i = 0; i < RK_HASHSIZE; i++)
  98     {
  99       target |= (Hashseq) NTN;
 100       if ((hashprobe & target) == target)
 101         return ((int) pos);
 102       target <<=4;
 103       pos++;
 104     }
 105
 106   return(-1);
 107 }
 108
 109
 110 #ifdef DEBUG                    /* Debugging aids */
 111
 112 static void
 113 writehash(Hashseq   hashseq)
 114 {
 115   int idx;
 116   int sym;
 117
 118   if (hashseq/16)
 119     writehash(hashseq/16);
 120
 121   sym = (int) (hashseq % 16);
 122   if (sym == 0)
 123     putchar('-');
 124   else
 125     {
 126       for (idx = 0; sym != iupac[idx].code && idx < IUPACSYMNUM; idx++);
 127       if (idx > IUPACSYMNUM)
 128         printf("(%d)", sym);
 129       else
 130         putchar(iupac[idx].sym);
 131     }
 132 }
 133
 134 #endif