src/fastq-sw.c

   1 /*
   2  * This file is part of fastq-tools.
   3  *
   4  * Copyright (c) 2011 by Daniel C. Jones <dcjones@cs.washington.edu>
   5  *
   6  * This is a very simple, 'fast enough' implementation of the Smith-Waterman
   7  * algorithm specifically for short nucleotide sequences, working in O(mn) time
   8  * and O(m) space, implemented according to the original Gotoh paper and
   9  * Phil Green's implementation in cross_match.
  10  *
  11  * There is no fancy bit packing or vectorization, but such features would offer
  12  * diminishing returns when aligning short sequences such as high throughput
  13  * sequencing data. For example the Farrar SSE2 algorithm might be a tiny bit
  14  * faster, but would diminish portability.
  15  *
  16  */
  17
  18
  19 #include "fastq-sw.h"
  20 #include "fastq-common.h"
  21 #include <string.h>
  22 #include <stdlib.h>
  23
  24
  25 static const int sw_default_d[25] =
  26     /* A   C   G   T   N */
  27     {  1, -2, -2, -2 , 0,
  28       -2,  1, -2, -2,  0,
  29       -2, -2,  1, -2,  0,
  30       -2, -2, -2,  1,  0,
  31        0,  0,  0,  0,  0  };
  32
  33 static const int sw_mat50_d[25] =
  34     {  2, -2,  0, -2 , 0,
  35       -2,  2, -2,  0,  0,
  36        0, -2,  2, -2,  0,
  37       -2,  0, -2,  2,  0,
  38        0,  0,  0,  0,  0  };
  39
  40 static const int sw_mat70_d[25] =
  41     {  2, -2, -1, -2 , 0,
  42       -2,  2, -2, -1,  0,
  43       -1, -2,  2, -2,  0,
  44       -2, -1, -2,  2,  0,
  45        0,  0,  0,  0,  0  };
  46
  47
  48 static inline int imax(int a, int b)
  49 {
  50     return a > b ? a : b;
  51 }
  52
  53 static inline int imax4(int a, int b, int c, int d)
  54 {
  55     return imax(imax(a,b), imax(c,d));
  56 }
  57
  58
  59
  60 void fastq_sw_conv_seq(unsigned char* seq, int n)
  61 {
  62     while (*seq || n--) {
  63         switch (*seq) {
  64             case 'A' :
  65             case 'a':
  66             case 'U':
  67             case 'u':
  68                 *seq = 0;
  69                 break;
  70
  71             case 'C':
  72             case 'c':
  73                 *seq = 1;
  74                 break;
  75
  76             case 'G':
  77             case 'g':
  78                 *seq = 2;
  79                 break;
  80
  81             case 'T':
  82             case 't':
  83                 *seq = 3;
  84                 break;
  85
  86             case 'N':
  87             case 'n':
  88             default:
  89                 *seq = 4;
  90         }
  91
  92         seq++;
  93     }
  94 }
  95
  96
  97 sw_t* fastq_alloc_sw(const unsigned char* subject, int size)
  98 {
  99     sw_t* sw = malloc_or_die(sizeof(sw_t));
 100     memcpy(sw->subject, subject, size);
 101
 102     /* default cost matrix */
 103     memcpy(sw->d, sw_default_d, 25 * sizeof(int));
 104
 105     /* default gap costs */
 106     sw->gap_open   = -4;
 107     sw->gap_extend = -3;
 108
 109     sw->work0 = malloc_or_die(size * sizeof(int));
 110     sw->work1 = malloc_or_die(size * sizeof(int));
 111     sw->size   = size;
 112
 113     return sw;
 114 }
 115
 116
 117 void fastq_free_sw(sw_t* sw)
 118 {
 119     free(sw->subject);
 120     free(sw->work0);
 121     free(sw->work1);
 122     free(sw);
 123 }
 124
 125
 126
 127 int fastq_sw(sw_t* sw, const unsigned char* x, int n)
 128 {
 129     /*const int max_score = 65535;*/
 130
 131     /* conveniance */
 132     int*      maxstu = sw->work0;
 133     int*           t = sw->work1;
 134     int            m = sw->size;
 135     const int*     d = sw->d;
 136     int       gap_op = sw->gap_open;
 137     int       gap_ex = sw->gap_extend;
 138     unsigned char* y = sw->subject;
 139
 140     int i, j;
 141
 142     /* zero maxstu row */
 143     memset(maxstu, 0, m * sizeof(int));
 144
 145     /* initialize t row to a negative value to prohibit
 146      * leading with gap extensions */
 147     for (j = 0; j < m; j++) t[j] = -1;
 148
 149     int u, s;
 150     int maxstu0;
 151
 152     for (i = 0; i < n; i++) {
 153
 154         /* special case for column 0 */
 155         t[j] = imax(maxstu[j] + gap_op, t[j] + gap_ex);
 156         u    = gap_op;
 157         maxstu[0] = imax4(0, d[5 * y[j] + x[j]], t[j], u);
 158         maxstu0 = 0;
 159
 160
 161         for (j = 1; j < m; j++) {
 162             t[j] = imax(maxstu[j] + gap_op, t[j] + gap_ex);
 163             u    = imax(maxstu[j-1] + gap_op, u + gap_ex);
 164             s    = maxstu0 + d[5 * y[j] + x[i]];
 165             maxstu0 = maxstu[j];
 166
 167             maxstu[j] = imax4(0, s, t[j], u);
 168         }
 169     }
 170
 171     /* return maximum from maxstu */
 172     maxstu0 = 0;
 173     for (j = 0; j <m; j++) maxstu0 = imax(maxstu0, maxstu[j]);
 174
 175     return maxstu0;
 176 }
 177
 178
 179