src/swsse2/swsse2.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 #include <string.h>\r
  17 \r
  18 #include <time.h>\r
  19 #include <sys/timeb.h>\r
  20 \r
  21 #include "swsse2.h"\r
  22 #include "matrix.h"\r
  23 #include "fastalib.h"\r
  24 \r
  25 #include "swscalar.h"\r
  26 #include "swwozniak.h"\r
  27 #ifdef WITH_ROGNES\r
  28 #include "swrognes.h"\r
  29 #endif\r
  30 #include "swstriped.h"\r
  31 \r
  32 typedef enum { SCALAR, \r
  33                WOZNIAK, \r
  34 #ifdef WITH_ROGNES\r
  35                ROGNES, \r
  36 #endif\r
  37                STRIPED \r
  38 } SW_TYPES;\r
  39 \r
  40 const char *SW_IMPLEMENATION[] = {\r
  41     "Non-optimized",\r
  42     "Wozniak",\r
  43 #ifdef WITH_ROGNES\r
  44     "Rognes",\r
  45 #endif\r
  46     "Striped",\r
  47 };\r
  48 \r
  49 typedef struct {\r
  50     SW_DATA *(*init) (unsigned char   *querySeq,\r
  51                       int              queryLength,\r
  52                       signed char     *matrix);\r
  53     void     (*scan) (unsigned char   *querySeq,\r
  54                       int              queryLength,\r
  55                       FASTA_LIB       *dbLib,\r
  56                       void            *swData,\r
  57                       SEARCH_OPTIONS  *options,\r
  58                       SCORE_LIST      *scores);\r
  59     void     (*done) (SW_DATA         *pSwData);\r
  60 } SW_FUNCT_DEFS;\r
  61 \r
  62 SW_FUNCT_DEFS swFuncs[] = {\r
  63     {\r
  64         swScalarInit,\r
  65         swScalarScan,\r
  66         swScalarComplete,\r
  67     },\r
  68     {\r
  69         swWozniakInit,\r
  70         swWozniakScan,\r
  71         swWozniakComplete,\r
  72     },\r
  73 #ifdef WITH_ROGNES\r
  74     {\r
  75         swRognesInit,\r
  76         swRognesScan,\r
  77         swRognesComplete,\r
  78     },\r
  79 #endif\r
  80     {\r
  81         swStripedInit,\r
  82         swStripedScan,\r
  83         swStripedComplete,\r
  84     },\r
  85 };\r
  86 \r
  87 const char AMINO_ACIDS[ALPHA_SIZE] = {\r
  88     'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', \r
  89     'I', 'K', 'L', 'M', 'N', 'P', 'Q', 'R', \r
  90     'S', 'T', 'V', 'W', 'X', 'Y', 'Z'\r
  91 };\r
  92 \r
  93 const int AMINO_ACID_VALUE[256] = {\r
  94     -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,\r
  95     -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,\r
  96     -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,\r
  97     -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,\r
  98     -1,  0,  1,  2,  3,  4,  5,  6,  7,  8, -1,  9, 10, 11, 12, -1,\r
  99     13, 14, 15, 16, 17, -1, 18, 19, 20, 21, 22, -1, -1, -1, -1, -1,\r
 100     -1,  0,  1,  2,  3,  4,  5,  6,  7,  8, -1,  9, 10, 11, 12, -1,\r
 101     13, 14, 15, 16, 17, -1, 18, 19, 20, 21, 22, -1, -1, -1, -1, -1,\r
 102     -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,\r
 103     -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,\r
 104     -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,\r
 105     -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,\r
 106     -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,\r
 107     -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,\r
 108     -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,\r
 109     -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1\r
 110 };\r
 111 \r
 112 SCORE_LIST *initList (int count);\r
 113 void freeList (SCORE_LIST *list);\r
 114 \r
 115 void printResults (SCORE_LIST *list);\r
 116 \r
 117 void printUsage (void)\r
 118 {\r
 119     printf ("Usage: swsse2 [-h] [-(n|w|r|s)] [-i num] [-e num] [-t num] "\r
 120         "[-c num] matrix query db\n");\r
 121     printf ("    -h       : this help message\n");\r
 122     printf ("    -n       : run a non-vectorized Smith-Waterman search\n");\r
 123     printf ("    -w       : run a vectorized Wozniak search\n");\r
 124 #ifdef WITH_ROGNES\r
 125     printf ("    -r       : run a vectorized Rognes search\n");\r
 126 #else\r
 127     printf ("    -r       : run a vectorized Rognes search (NOT SUPPORTED)\n");\r
 128 #endif\r
 129     printf ("    -s       : run a vectorized striped search (default)\n");\r
 130     printf ("    -i num   : gap init penalty (default -10)\n");\r
 131     printf ("    -e num   : gap extension penalty (default -2)\n");\r
 132     printf ("    -t num   : minimum score threshold (default 20)\n");\r
 133     printf ("    -c num   : number of scores to be displayed (default 250)\n");\r
 134     printf ("    matrix   : scoring matrix file\n");\r
 135     printf ("    query    : query sequence file (fasta format)\n");\r
 136     printf ("    db       : sequence database file (fasta format)\n");\r
 137 }\r
 138 \r
 139 #if 0\r
 140 int main (int argc, char **argv)\r
 141 {\r
 142     int i;\r
 143 \r
 144     int penalty;\r
 145 \r
 146     int rptCount = 250;\r
 147 \r
 148     SW_TYPES swType = STRIPED;\r
 149 \r
 150     char *dbFile = NULL;\r
 151     char *queryFile = NULL;\r
 152     char *matrixFile = NULL;\r
 153 \r
 154     signed char *matrix;\r
 155 \r
 156     unsigned char *querySeq;\r
 157     int queryLen;\r
 158 \r
 159     SCORE_LIST *list;\r
 160 \r
 161     FASTA_LIB *queryLib;\r
 162     FASTA_LIB *dbLib;\r
 163 \r
 164     void *swData;\r
 165 \r
 166     struct _timeb startTime;\r
 167     struct _timeb endTime;\r
 168 \r
 169     double startTimeSec;\r
 170     double endTimeSec;\r
 171 \r
 172     SEARCH_OPTIONS options;\r
 173 \r
 174     if (argc < 4) {\r
 175         printUsage ();\r
 176         exit (-1);\r
 177     }\r
 178 \r
 179     /* set the default options */\r
 180     options.gapInit = -10;\r
 181     options.gapExt = -2;\r
 182     options.threshold = 20;\r
 183 \r
 184     i = 1;\r
 185     while (i < argc) {\r
 186         if (i + 3 == argc) {\r
 187             /* should be matrix file name */\r
 188             matrixFile = argv[i];\r
 189 \r
 190         } else if (i + 2 == argc) {\r
 191             /* should be query file name */\r
 192             queryFile = argv[i];\r
 193 \r
 194         } else if (i + 1 == argc) {\r
 195             /* should be matrix file name */\r
 196             dbFile = argv[i];\r
 197 \r
 198         } else {\r
 199             /* process arguements */\r
 200             switch (argv[i][1]) {\r
 201                case 'h':\r
 202                    printUsage ();\r
 203                    break;\r
 204                case 'n':\r
 205                    swType = SCALAR;\r
 206                    break;\r
 207                case 'r':\r
 208 #ifdef WITH_ROGNES\r
 209                    swType = ROGNES;\r
 210 #else\r
 211                    fprintf (stderr, "The ROGNES implementation is not supported\n");\r
 212                    exit (-1);\r
 213 #endif\r
 214                    break;\r
 215                case 'w':\r
 216                    swType = WOZNIAK;\r
 217                    break;\r
 218                case 's':\r
 219                    swType = STRIPED;\r
 220                    break;\r
 221                case 'i':\r
 222                    penalty = atoi (argv[++i]);\r
 223                    if (penalty > 0 || penalty < -128) {\r
 224                         fprintf (stderr, "Invalid gap init %d\n", penalty);\r
 225                         fprintf (stderr, "Valid range is 0 - -128\n");\r
 226                         exit (-1);\r
 227                    }\r
 228                    options.gapInit = (unsigned short) penalty;\r
 229                    break;\r
 230                case 'e':\r
 231                    penalty = atoi (argv[++i]);\r
 232                    if (penalty > 0 || penalty < -128) {\r
 233                         fprintf (stderr, "Invalid gap extension %d\n", penalty);\r
 234                         fprintf (stderr, "Valid range is 0 - -128\n");\r
 235                         exit (-1);\r
 236                    }\r
 237                    options.gapExt = (unsigned short) penalty;\r
 238                    break;\r
 239                case 't':\r
 240                    options.threshold = atoi (argv[++i]);\r
 241                    break;\r
 242                case 'c':\r
 243                    rptCount = atoi (argv[++i]);\r
 244                    if (rptCount < 10) {\r
 245                        rptCount = 10;\r
 246                    }\r
 247                    break;\r
 248                default:\r
 249                    fprintf (stderr, "Invalid option %s\n", argv[i]);\r
 250                    printUsage ();\r
 251                    exit (-1);\r
 252             }\r
 253         }\r
 254         ++i;\r
 255     }\r
 256 \r
 257     list = initList (rptCount);\r
 258 \r
 259     matrix = readMatrix (matrixFile);\r
 260     if (matrix == NULL) {\r
 261         fprintf (stderr, "Error reading matrix\n");\r
 262         exit (-1);\r
 263     }\r
 264 \r
 265     dbLib = openLib (dbFile, swType == WOZNIAK);\r
 266     queryLib = openLib (queryFile, 0);\r
 267 \r
 268     querySeq = nextSeq (queryLib, &queryLen);\r
 269     if (queryLen == 0) {\r
 270         fprintf (stderr, "Empty query sequence\n");\r
 271         exit (-1);\r
 272     }\r
 273 \r
 274     printf ("%s vs %s\n", queryFile, dbFile);\r
 275     printf ("Matrix: %s, Init: %d, Ext: %d\n\n", \r
 276             matrixFile, options.gapInit, options.gapExt);\r
 277 \r
 278     _ftime(&startTime);\r
 279 \r
 280     swData = (swFuncs[swType].init) (querySeq, queryLen, matrix);\r
 281 \r
 282     (swFuncs[swType].scan) (querySeq, queryLen, dbLib, swData, &options, list);\r
 283 \r
 284     (swFuncs[swType].done) (swData);\r
 285 \r
 286     _ftime(&endTime);\r
 287 \r
 288     printResults (list);\r
 289 \r
 290     printf ("\n");\r
 291     printf ("%d residues in query string\n", queryLen);\r
 292     printf ("%d residues in %d library sequences\n", \r
 293             dbLib->residues, dbLib->sequences);\r
 294 \r
 295     startTimeSec = startTime.time + startTime.millitm / 1000.0;\r
 296     endTimeSec = endTime.time + endTime.millitm / 1000.0;\r
 297     printf ("Scan time: %6.3f (%s implementation)\n", \r
 298             endTimeSec - startTimeSec,\r
 299             SW_IMPLEMENATION[swType]);\r
 300   \r
 301     closeLib (queryLib);\r
 302     closeLib (dbLib);\r
 303 \r
 304     free (matrix);\r
 305 \r
 306     freeList (list);\r
 307 }\r
 308 #endif\r
 309 \r
 310 SCORE_LIST *\r
 311 initList (int count)\r
 312 {\r
 313     int i;\r
 314 \r
 315     SCORE_LIST *hdr;\r
 316     SCORE_NODE *list;\r
 317     SCORE_NODE *prev;\r
 318 \r
 319     hdr = (SCORE_LIST *) malloc (sizeof (SCORE_LIST));\r
 320     if (hdr == NULL) {\r
 321         fprintf (stderr, "Cannot allocate storage for score header\n");\r
 322         exit (-1);\r
 323     }\r
 324 \r
 325     list = (SCORE_NODE *) calloc (count, sizeof (SCORE_NODE));\r
 326     if (list == NULL) {\r
 327         fprintf (stderr, "Cannot allocate storage for scores\n");\r
 328         exit (-1);\r
 329     }\r
 330 \r
 331     /* initialize the scores list */\r
 332     hdr->minScore = 0;\r
 333     hdr->first = NULL;\r
 334     hdr->last = NULL;\r
 335     hdr->free = list;\r
 336     hdr->buffer = list;\r
 337 \r
 338     prev = NULL;\r
 339     for (i = 0; i < count; ++i) {\r
 340         list[i].name[0] = '\0';\r
 341         list[i].score = 0;\r
 342 \r
 343         if (i == 0) {\r
 344             list[i].prev = NULL;\r
 345         } else {\r
 346             list[i].prev = &list[i-1];\r
 347         }\r
 348 \r
 349         if (i == count - 1) {\r
 350             list[i].next = NULL;\r
 351         } else {\r
 352             list[i].next = &list[i+1];\r
 353         }\r
 354     }\r
 355 \r
 356     return hdr;\r
 357 }\r
 358 \r
 359 void freeList (SCORE_LIST *list)\r
 360 {\r
 361     free (list->buffer);\r
 362     free (list);\r
 363 }\r
 364 \r
 365 int insertList (SCORE_LIST *list, int score, char *name)\r
 366 {\r
 367     SCORE_NODE *node;\r
 368     SCORE_NODE *ptr = list->first;\r
 369 \r
 370     if (list->free != NULL) {\r
 371         node = list->free;\r
 372         list->free = list->free->next;\r
 373     } else if (score > list->last->score) {\r
 374         node = list->last;\r
 375         list->last = node->prev;\r
 376         list->last->next = NULL;\r
 377     } else {\r
 378         /* should never happen */\r
 379         return list->minScore + 1;\r
 380     }\r
 381 \r
 382     strncpy (node->name, name, MAX_SCORE_NAME);\r
 383     node->name[MAX_SCORE_NAME - 1] = '\0';\r
 384     node->score = score;\r
 385 \r
 386     while (ptr && ptr->score >= score) {\r
 387         ptr = ptr->next;\r
 388     }\r
 389 \r
 390     if (list->first == NULL) {\r
 391         list->first = node;\r
 392         list->last = node;\r
 393         node->prev = NULL;\r
 394         node->next = NULL;\r
 395     } else if (ptr == NULL) {\r
 396         node->prev = list->last;\r
 397         node->next = NULL;\r
 398         node->prev->next  = node;\r
 399         list->last = node;\r
 400     } else {\r
 401         node->prev = ptr->prev;\r
 402         node->next = ptr;\r
 403 \r
 404         if (node->prev == NULL) {\r
 405             list->first = node;\r
 406         } else {\r
 407             node->prev->next = node;\r
 408         }\r
 409         ptr->prev = node;\r
 410     }\r
 411 \r
 412     if (list->free == NULL) {\r
 413         list->minScore = list->last->score + 1;\r
 414     }\r
 415 \r
 416     return list->minScore;\r
 417 }\r
 418 \r
 419 void printResults (SCORE_LIST *list)\r
 420 {\r
 421     SCORE_NODE *ptr = list->first;\r
 422 \r
 423     printf ("Score  Description\n");\r
 424 \r
 425     while (ptr) {\r
 426         printf ("%5d  %s\n", ptr->score, ptr->name);\r
 427         ptr = ptr->next;\r
 428     }\r
 429 }\r
 430 \r