src/swsse2/swscalar.c

   1 /******************************************************************\r
   2   Copyright 2006 by Michael Farrar.  All rights reserved.\r
   3   This program may not be sold or incorporated into a commercial product,\r
   4   in whole or in part, without written consent of Michael Farrar.  For \r
   5   further information regarding permission for use or reproduction, please \r
   6   contact: Michael Farrar at farrar.michael@gmail.com.\r
   7 *******************************************************************/\r
   8 \r
   9 /*\r
  10   Written by Michael Farrar, 2006.\r
  11   Please send bug reports and/or suggestions to farrar.michael@gmail.com.\r
  12 */\r
  13 \r
  14 #include <stdio.h>\r
  15 #include <stdlib.h>\r
  16 \r
  17 #include "swsse2.h"\r
  18 #include "swscalar.h"\r
  19 \r
  20 typedef struct {\r
  21     int  *pMatrix;\r
  22     int  *pH;\r
  23     int  *pE;\r
  24     int  *pData;\r
  25 } SwScalarData;\r
  26 \r
  27 int\r
  28 swScalar (unsigned char   *querySeq,\r
  29           int              queryLength,\r
  30           unsigned char   *dbSeq,\r
  31           int              dbLength,\r
  32           unsigned short   gapOpen,\r
  33           unsigned short   gapExtend,\r
  34           int             *pMatrix,\r
  35           int             *pH,\r
  36           int             *pE);\r
  37 \r
  38 void *\r
  39 swScalarInit(unsigned char   *querySeq,\r
  40              int              queryLength,\r
  41              signed char     *matrix)\r
  42 {\r
  43     int i, j;\r
  44     int nCount;\r
  45 \r
  46     int *pi;\r
  47 \r
  48     SwScalarData *pSwData;\r
  49  \r
  50     /* remove unreferenced warnings */\r
  51     querySeq;\r
  52 \r
  53     pSwData = (SwScalarData *) malloc (sizeof (SwScalarData));\r
  54     if (!pSwData) {\r
  55         fprintf (stderr, "Unable to allocate memory for SW data\n");\r
  56         exit (-1);\r
  57     }\r
  58 \r
  59     nCount = ALPHA_SIZE * ALPHA_SIZE +        /* matrix */\r
  60              (queryLength * 3);               /* vH1, vH2 and vE */\r
  61 \r
  62     pSwData->pData = (int *) calloc (nCount, sizeof (int));\r
  63     if (!pSwData->pData) {\r
  64         fprintf (stderr, "Unable to allocate memory for SW data buffers\n");\r
  65         exit (-1);\r
  66     }\r
  67 \r
  68     pSwData->pMatrix = pSwData->pData;\r
  69 \r
  70     pSwData->pH = pSwData->pMatrix + ALPHA_SIZE * ALPHA_SIZE;\r
  71     pSwData->pE  = pSwData->pH + queryLength;\r
  72 \r
  73     /* Conver the scoring matrix to ints */\r
  74     pi = pSwData->pMatrix;\r
  75     for (i = 0; i < ALPHA_SIZE; ++i) {\r
  76         for (j = 0; j < ALPHA_SIZE; ++j) {\r
  77             *pi++ = (int) *matrix++;\r
  78         }\r
  79     }\r
  80 \r
  81     return pSwData;\r
  82 }\r
  83 \r
  84 void swScalarScan (unsigned char   *querySeq,\r
  85                    int              queryLength,\r
  86                    FASTA_LIB       *dbLib,\r
  87                    void            *swData,\r
  88                    SEARCH_OPTIONS  *options,\r
  89                    SCORE_LIST      *scores)\r
  90 {\r
  91     int score;\r
  92 \r
  93     int threshold = options->threshold;\r
  94 \r
  95     unsigned char *dbSeq;\r
  96     int dbLen;\r
  97 \r
  98     int gapInit = -(options->gapInit + options->gapExt);\r
  99     int gapExt  = -options->gapExt;\r
 100 \r
 101     SwScalarData *scalarData = (SwScalarData *) swData;\r
 102 \r
 103     dbSeq = nextSeq (dbLib, &dbLen);\r
 104     while (dbLen > 0) {\r
 105 \r
 106         score = swScalar (querySeq, queryLength, \r
 107                           dbSeq, dbLen, \r
 108                           gapInit, gapExt, \r
 109                           scalarData->pMatrix,\r
 110                           scalarData->pH,\r
 111                           scalarData->pE);\r
 112 \r
 113         if (score >= threshold) {\r
 114             int minScore = insertList (scores, score, seqName (dbLib));\r
 115             if (minScore >= threshold) {\r
 116                 threshold = minScore;\r
 117             }\r
 118         }\r
 119 \r
 120         dbSeq = nextSeq (dbLib, &dbLen);\r
 121     }\r
 122 }\r
 123 \r
 124 void\r
 125 swScalarComplete(void *pSwData)\r
 126 {\r
 127     SwScalarData *pScalarData = (SwScalarData *) pSwData;\r
 128 \r
 129     free (pScalarData->pData);\r
 130     free (pScalarData);\r
 131 }\r
 132 \r
 133 int\r
 134 swScalar(unsigned char   *querySeq,\r
 135          int              queryLength,\r
 136          unsigned char   *dbSeq,\r
 137          int              dbLength,\r
 138          unsigned short   gapOpen,\r
 139          unsigned short   gapExtend,\r
 140          int             *pMatrix,\r
 141          int             *pH,\r
 142          int             *pE)\r
 143 {\r
 144     int     i, j;\r
 145 \r
 146     int E, F, H;\r
 147     int prevH;\r
 148 \r
 149     int maxScore = 0;\r
 150 \r
 151     int *pScore;\r
 152 \r
 153     /* Zero out the storage vector */\r
 154     for (i = 0; i < queryLength; i++)\r
 155     {\r
 156         *(pE + i) = 0;\r
 157         *(pH + i) = 0;\r
 158     }\r
 159 \r
 160     for (i = 0; i < dbLength; ++i)\r
 161     {\r
 162         /* fetch first data asap. */\r
 163         pScore = pMatrix + dbSeq[i] * ALPHA_SIZE;\r
 164 \r
 165         /* zero out F. */\r
 166         F = 0;\r
 167 \r
 168         /* load the next h value */\r
 169         prevH = 0;\r
 170 \r
 171         for (j = 0; j < queryLength; j++)\r
 172         {\r
 173             /* load E values */\r
 174             E = *(pE + j);\r
 175 \r
 176             /* add score to H */\r
 177             H = prevH + *(pScore + querySeq[j]);\r
 178             if (H < 0) {\r
 179                 H = 0;\r
 180             }\r
 181 \r
 182             /* Update highest score encountered this far */\r
 183             if (H > maxScore) {\r
 184                 maxScore = H;\r
 185             }\r
 186 \r
 187             /* get max from H, E and F */\r
 188             if (E > H) {\r
 189                 H = E;\r
 190             }\r
 191             if (F > H) {\r
 192                 H = F;\r
 193             }\r
 194 \r
 195             /* save H values */\r
 196             prevH = *(pH + j);\r
 197             *(pH + j) = H;\r
 198 \r
 199             H = H - gapOpen;\r
 200 \r
 201             /* update E value */\r
 202             E = E - gapExtend;\r
 203             if (H > E) {\r
 204                 E = H;\r
 205             }\r
 206             if (E < 0) {\r
 207                 E = 0;\r
 208             }\r
 209 \r
 210             /* update vF value */\r
 211             F = F - gapExtend;\r
 212             if (H > F) {\r
 213                 F = H;\r
 214             }\r
 215             if (F < 0) {\r
 216                 F = 0;\r
 217             }\r
 218 \r
 219             /* save vE values */\r
 220             *(pE + j) = E;\r
 221         }\r
 222     }\r
 223 \r
 224     /* return largest score */\r
 225     return maxScore;\r
 226 }\r