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         n--;
  94     }
  95 }
  96
  97
  98 sw_t* fastq_alloc_sw(const unsigned char* subject, int size)
  99 {
 100     sw_t* sw = malloc_or_die(sizeof(sw_t));
 101
 102     sw->subject = malloc_or_die(size);
 103     memcpy(sw->subject, subject, size);
 104
 105     /* default cost matrix */
 106     memcpy(sw->d, sw_default_d, 25 * sizeof(int));
 107
 108     /* default gap costs */
 109     sw->gap_open   = -4;
 110     sw->gap_extend = -3;
 111
 112     sw->work0 = malloc_or_die(size * sizeof(int));
 113     sw->work1 = malloc_or_die(size * sizeof(int));
 114     sw->size   = size;
 115
 116     return sw;
 117 }
 118
 119
 120 void fastq_free_sw(sw_t* sw)
 121 {
 122     free(sw->subject);
 123     free(sw->work0);
 124     free(sw->work1);
 125     free(sw);
 126 }
 127
 128
 129
 130 int fastq_sw(sw_t* sw, const unsigned char* x, int n)
 131 {
 132     /* conveniance */
 133     int*      maxstu = sw->work0;
 134     int*           t = sw->work1;
 135     int            m = sw->size;
 136     const int*     d = sw->d;
 137     int       gap_op = sw->gap_open;
 138     int       gap_ex = sw->gap_extend;
 139     unsigned char* y = sw->subject;
 140
 141     int i, j;
 142
 143     int score = 0;
 144
 145     /* zero maxstu row */
 146     memset(maxstu, 0, m * sizeof(int));
 147
 148     /* initialize t row to a negative value to prohibit
 149      * leading with gap extensions */
 150     for (j = 0; j < m; j++) t[j] = -1;
 151
 152     int u, s;
 153     int maxstu0;
 154
 155     for (i = 0; i < n; i++) {
 156
 157         /* special case for column 0 */
 158         t[0] = imax(maxstu[0] + gap_op, t[0] + gap_ex);
 159         u    = gap_op;
 160         maxstu[0] = imax4(0, d[5 * y[0] + x[0]], t[0], u);
 161         maxstu0 = 0;
 162         score = imax(score, maxstu[0]);
 163
 164
 165         for (j = 1; j < m; j++) {
 166             t[j] = imax(maxstu[j] + gap_op, t[j] + gap_ex);
 167             u    = imax(maxstu[j-1] + gap_op, u + gap_ex);
 168             s    = maxstu0 + d[5 * y[j] + x[i]];
 169             maxstu0 = maxstu[j];
 170
 171             maxstu[j] = imax4(0, s, t[j], u);
 172             score = imax(score, maxstu[j]);
 173         }
 174     }
 175
 176     return score;
 177 }
 178
 179
 180