2 * This file is part of fastq-tools.
4 * Copyright (c) 2011 by Daniel C. Jones <dcjones@cs.washington.edu>
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.
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.
25 static const int sw_default_d[25] =
33 static const int sw_mat50_d[25] =
40 static const int sw_mat70_d[25] =
48 static inline int imax(int a, int b)
53 static inline int imax4(int a, int b, int c, int d)
55 return imax(imax(a,b), imax(c,d));
60 void fastq_sw_conv_seq(unsigned char* seq, int n)
98 sw_t* fastq_alloc_sw(const unsigned char* subject, int size)
100 sw_t* sw = malloc_or_die(sizeof(sw_t));
102 sw->subject = malloc_or_die(size);
103 memcpy(sw->subject, subject, size);
105 /* default cost matrix */
106 memcpy(sw->d, sw_default_d, 25 * sizeof(int));
108 /* default gap costs */
112 sw->work0 = malloc_or_die(size * sizeof(int));
113 sw->work1 = malloc_or_die(size * sizeof(int));
120 void fastq_free_sw(sw_t* sw)
130 int fastq_sw(sw_t* sw, const unsigned char* x, int n)
133 int* maxstu = sw->work0;
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;
145 /* zero maxstu row */
146 memset(maxstu, 0, m * sizeof(int));
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;
155 for (i = 0; i < n; i++) {
157 /* special case for column 0 */
158 t[0] = imax(maxstu[0] + gap_op, t[0] + gap_ex);
160 maxstu[0] = imax4(0, d[5 * y[0] + x[0]], t[0], u);
162 score = imax(score, maxstu[0]);
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]];
171 maxstu[j] = imax4(0, s, t[j], u);
172 score = imax(score, maxstu[j]);