Branch data Line data Source code
1 : : // $Id: SmithWaterman.cc 6219 2008-10-01 05:39:07Z vern $
2 : : //
3 : : // See the file "COPYING" in the main distribution directory for copyright.
4 : :
5 : : #include "config.h"
6 : :
7 : : #include <algorithm>
8 : : #include <ctype.h>
9 : :
10 : : #include "SmithWaterman.h"
11 : : #include "Var.h"
12 : : #include "util.h"
13 : :
14 : 0 : BroSubstring::BroSubstring(const BroSubstring& bst)
15 : 0 : : BroString((const BroString&) bst), _new(bst._new)
16 : : {
17 [ # # ][ # # ]: 0 : for ( BSSAlignVecCIt it = bst._aligns.begin(); it != bst._aligns.end(); ++it )
18 : 0 : _aligns.push_back(*it);
19 : 0 : }
20 : :
21 : 0 : const BroSubstring& BroSubstring::operator=(const BroSubstring& bst)
22 : : {
23 : 0 : BroString::operator=(bst);
24 : :
25 : 0 : _aligns.clear();
26 : :
27 [ # # ]: 0 : for ( BSSAlignVecCIt it = bst._aligns.begin(); it != bst._aligns.end(); ++it )
28 : 0 : _aligns.push_back(*it);
29 : :
30 : 0 : _new = bst._new;
31 : :
32 : 0 : return *this;
33 : : }
34 : :
35 : 0 : void BroSubstring::AddAlignment(const BroString* str, int index)
36 : : {
37 : 0 : _aligns.push_back(BSSAlign(str, index));
38 : 0 : }
39 : :
40 : 0 : bool BroSubstring::DoesCover(const BroSubstring* bst) const
41 : : {
42 [ # # ]: 0 : if ( _aligns.size() != bst->_aligns.size() )
43 : 0 : return false;
44 : :
45 : 0 : BSSAlignVecCIt it_bst = bst->_aligns.begin();
46 : :
47 [ # # ]: 0 : for ( BSSAlignVecCIt it = _aligns.begin(); it != _aligns.end(); ++it, ++it_bst )
48 : : {
49 : 0 : const BSSAlign& a = *it;
50 : 0 : const BSSAlign& a_bst = *it_bst;
51 : :
52 [ # # # # ]: 0 : if (a.index > a_bst.index || a.index + Len() < a_bst.index + bst->Len())
[ # # ]
53 : 0 : return false;
54 : : }
55 : :
56 : 0 : return true;
57 : : }
58 : :
59 : 0 : VectorVal* BroSubstring::VecToPolicy(Vec* vec)
60 : : {
61 : : RecordType* sw_substring_type =
62 : 0 : internal_type("sw_substring")->AsRecordType();
63 [ # # ]: 0 : if ( ! sw_substring_type )
64 : 0 : return 0;
65 : :
66 : : RecordType* sw_align_type =
67 : 0 : internal_type("sw_align")->AsRecordType();
68 [ # # ]: 0 : if ( ! sw_align_type )
69 : 0 : return 0;
70 : :
71 : : VectorType* sw_align_vec_type =
72 : 0 : internal_type("sw_align_vec")->AsVectorType();
73 [ # # ]: 0 : if ( ! sw_align_vec_type )
74 : 0 : return 0;
75 : :
76 : : VectorVal* result =
77 : 0 : new VectorVal(internal_type("sw_substring_vec")->AsVectorType());
78 [ # # ]: 0 : if ( ! result )
79 : 0 : return 0;
80 : :
81 [ # # ]: 0 : if ( vec )
82 : : {
83 [ # # ]: 0 : for ( unsigned int i = 0; i < vec->size(); ++i )
84 : : {
85 : 0 : BroSubstring* bst = (*vec)[i];
86 : :
87 : 0 : RecordVal* st_val = new RecordVal(sw_substring_type);
88 : 0 : st_val->Assign(0, new StringVal(new BroString(*bst)));
89 : :
90 : 0 : VectorVal* aligns = new VectorVal(sw_align_vec_type);
91 : :
92 [ # # ]: 0 : for ( unsigned int j = 0; j < bst->GetNumAlignments(); ++j )
93 : : {
94 : 0 : const BSSAlign& align = (bst->GetAlignments())[j];
95 : :
96 : 0 : RecordVal* align_val = new RecordVal(sw_align_type);
97 : 0 : align_val->Assign(0, new StringVal(new BroString(*align.string)));
98 : 0 : align_val->Assign(1, new Val(align.index, TYPE_COUNT));
99 : :
100 : 0 : aligns->Assign(j+1, align_val, 0);
101 : : }
102 : :
103 : 0 : st_val->Assign(1, aligns);
104 : 0 : st_val->Assign(2, new Val(bst->IsNewAlignment(), TYPE_BOOL));
105 : 0 : result->Assign(i+1, st_val, 0);
106 : : }
107 : : }
108 : :
109 : 0 : return result;
110 : : }
111 : :
112 : 0 : BroSubstring::Vec* BroSubstring::VecFromPolicy(VectorVal* vec)
113 : : {
114 : 0 : Vec* result = new Vec();
115 : :
116 : : // VectorVals start at index 1!
117 [ # # ]: 0 : for ( unsigned int i = 1; i <= vec->Size(); ++i )
118 : : {
119 : 0 : Val* v = vec->Lookup(i); // get the RecordVal
120 [ # # ]: 0 : if ( ! v )
121 : 0 : continue;
122 : :
123 : 0 : const BroString* str = v->AsRecordVal()->Lookup(0)->AsString();
124 : 0 : BroSubstring* substr = new BroSubstring(*str);
125 : :
126 : 0 : const VectorVal* aligns = v->AsRecordVal()->Lookup(1)->AsVectorVal();
127 [ # # ]: 0 : for ( unsigned int j = 1; j <= aligns->Size(); ++j )
128 : : {
129 : 0 : const RecordVal* align = aligns->AsVectorVal()->Lookup(j)->AsRecordVal();
130 : 0 : const BroString* str = align->Lookup(0)->AsString();
131 : 0 : int index = align->Lookup(1)->AsCount();
132 : 0 : substr->AddAlignment(str, index);
133 : : }
134 : :
135 : 0 : bool new_alignment = v->AsRecordVal()->Lookup(2)->AsBool();
136 : 0 : substr->MarkNewAlignment(new_alignment);
137 : :
138 : 0 : result->push_back(substr);
139 : : }
140 : :
141 : 0 : return result;
142 : : }
143 : :
144 : 0 : char* BroSubstring::VecToString(Vec* vec)
145 : : {
146 : 0 : string result("[");
147 : :
148 [ # # ]: 0 : for ( BroSubstring::VecIt it = vec->begin(); it != vec->end(); ++it )
149 : : {
150 : 0 : result += (*it)->CheckString();
151 : 0 : result += ",";
152 : : }
153 : :
154 : 0 : result += "]";
155 : 0 : return strdup(result.c_str());
156 : : }
157 : :
158 : 0 : BroString::IdxVec* BroSubstring::GetOffsetsVec(const Vec* vec, unsigned int index)
159 : : {
160 : 0 : BroString::IdxVec* result = new BroString::IdxVec();
161 : :
162 [ # # ]: 0 : for ( VecCIt it = vec->begin(); it != vec->end(); ++it )
163 : : {
164 : : int start, end;
165 : 0 : const BroSubstring* bst = (*it);
166 : :
167 [ # # ]: 0 : if ( bst->_aligns.size() <= index )
168 : 0 : continue;
169 : :
170 : 0 : const BSSAlign& align = bst->_aligns[index];
171 : 0 : start = align.index;
172 : 0 : end = start + bst->Len();
173 : :
174 : 0 : result->push_back(start);
175 : 0 : result->push_back(end);
176 : : }
177 : :
178 : 0 : return result;
179 : : }
180 : :
181 : :
182 : : bool BroSubstringCmp::operator()(const BroSubstring* bst1,
183 : 0 : const BroSubstring* bst2) const
184 : : {
185 [ # # ][ # # ]: 0 : if ( _index >= bst1->GetNumAlignments() ||
[ # # ]
186 : : _index >= bst2->GetNumAlignments() )
187 : : {
188 : 0 : warn("BroSubstringCmp::operator(): invalid index for input strings.\n");
189 : 0 : return false;
190 : : }
191 : :
192 [ # # ]: 0 : if ( bst1->GetAlignments()[_index].index <=
193 : : bst2->GetAlignments()[_index].index )
194 : 0 : return true;
195 : :
196 : 0 : return false;
197 : : }
198 : :
199 : : // A node in Smith-Waterman's dynamic programming matrix. Each node
200 : : // contains the byte it represents in the case of a match, the score
201 : : // at this point, and a pointer to the previous cell. Previous means
202 : : // one up and left in case of a match, or a jump somewhere above and
203 : : // left in case of a gap.
204 : : //
205 : : struct SWNode {
206 : : // ID field for the cell, for debugging purposes.
207 : : int id;
208 : :
209 : : u_char swn_byte;
210 : : bool swn_byte_assigned;
211 : : bool swn_visited;
212 : :
213 : : // The score in this cell. The cell with the globally best score
214 : : // marks the end of the alignment.
215 : : int swn_score;
216 : :
217 : : // Pointer to previous match, walking back yields subsequence.
218 : : SWNode* swn_prev;
219 : : };
220 : :
221 : : // A matrix of Smith-Waterman nodes.
222 : : //
223 : : class SWNodeMatrix {
224 : : public:
225 : 0 : SWNodeMatrix(const BroString* s1, const BroString* s2)
226 : 0 : : _s1(s1), _s2(s2), _rows(s1->Len() + 1), _cols(s2->Len() + 1)
227 : : {
228 : 0 : _nodes = new SWNode[_cols * _rows];
229 : 0 : memset(_nodes, 0, sizeof(SWNode) * _cols * _rows);
230 : 0 : }
231 : :
232 [ # # ]: 0 : ~SWNodeMatrix() { delete [] _nodes; }
233 : :
234 : 0 : SWNode* operator()(int row, int col)
235 : : {
236 : : // Make sure access is in allowed range.
237 [ # # ][ # # ]: 0 : if ( row < 0 || row >= _rows )
238 : 0 : return 0;
239 [ # # ][ # # ]: 0 : if ( col < 0 || col >= _cols )
240 : 0 : return 0;
241 : :
242 : 0 : return &(_nodes[row * _cols + col]);
243 : : }
244 : :
245 : 0 : const BroString* GetRowsString() const { return _s1; }
246 : 0 : const BroString* GetColsString() const { return _s2; }
247 : :
248 : 0 : int GetHeight() const { return _rows; }
249 : 0 : int GetWidth() const { return _cols; }
250 : :
251 : : // Quick helper function that calculates the coordinates of a
252 : : // node in the matrix via pointer arithmetic.
253 : : //
254 : 0 : void GetNodeIndices(SWNode* node, int& row, int& col)
255 : : {
256 : 0 : SWNode* base = &_nodes[0];
257 : 0 : int offset = (node - base);
258 : 0 : col = (offset % _cols);
259 : 0 : row = (offset / _cols);
260 : 0 : }
261 : :
262 : : private:
263 : : const BroString* _s1;
264 : : const BroString* _s2;
265 : :
266 : : int _rows, _cols;
267 : : SWNode* _nodes;
268 : : };
269 : :
270 : : // Returns the common subsequence starting from a given node.
271 : : // @result: vector holding results on return.
272 : : // @matrix: SW matrix.
273 : : // @node: starting node.
274 : : // @params: SW parameters.
275 : : //
276 : : static void sw_collect_single(BroSubstring::Vec* result, SWNodeMatrix& matrix,
277 : 0 : SWNode* node, SWParams& params)
278 : : {
279 : 0 : string substring("");
280 : 0 : int row = 0, col = 0;
281 : :
282 [ # # ]: 0 : while ( node )
283 : : {
284 : : // printf("NODE: %i\n", node->id);
285 : 0 : node->swn_visited = true;
286 : :
287 : : // Once we hit a gap, terminate the string and prepend
288 : : // it to our result vector, IF it has at least the length
289 : : // requested through the params._min_toklen parameter.
290 : : //
291 [ # # ]: 0 : if ( node->swn_byte_assigned )
292 : : {
293 : 0 : matrix.GetNodeIndices(node, row, col);
294 : 0 : substring += node->swn_byte;
295 : : // printf("SUBSTRING: %s\n", substring.c_str());
296 : : }
297 : : else
298 : : {
299 : : // printf("GAP\n");
300 [ # # ]: 0 : if ( substring.size() >= params._min_toklen )
301 : : {
302 : 0 : reverse(substring.begin(), substring.end());
303 : 0 : BroSubstring* bst = new BroSubstring(substring);
304 : 0 : bst->AddAlignment(matrix.GetRowsString(), row-1);
305 : 0 : bst->AddAlignment(matrix.GetColsString(), col-1);
306 : 0 : result->push_back(bst);
307 : : }
308 : :
309 : 0 : substring = "";
310 : : }
311 : :
312 : 0 : node = node->swn_prev;
313 : : }
314 : :
315 : : // Anything left over now is the first string of an alignment and is
316 : : // manually added and marked as the beginning of a new alignment.
317 : : //
318 [ # # ]: 0 : if ( substring.size() > 0 )
319 : : {
320 : 0 : reverse(substring.begin(), substring.end());
321 : 0 : BroSubstring* bst = new BroSubstring(substring);
322 : 0 : bst->AddAlignment(matrix.GetRowsString(), row-1);
323 : 0 : bst->AddAlignment(matrix.GetColsString(), col-1);
324 : 0 : result->push_back(bst);
325 : : }
326 : :
327 [ # # ]: 0 : if ( result->size() > 0 )
328 : 0 : result->back()->MarkNewAlignment(true);
329 : 0 : }
330 : :
331 : : // Returns repeated common-subsequence alignments.
332 : : // @result: vector holding results on return.
333 : : // @matrix: SW matrix.
334 : : // @params: SW parameters.
335 : : //
336 : : // The approach taken is to essentially follow back from all starting points of
337 : : // common subsequences while tracking which nodes were visited earlier and which
338 : : // substrings are redundant (i.e., fully covered by a larger common substring).
339 : : //
340 : : static void sw_collect_multiple(BroSubstring::Vec* result,
341 : 0 : SWNodeMatrix& matrix, SWParams& params)
342 : : {
343 : 0 : vector<BroSubstring::Vec*> als;
344 : :
345 [ # # ]: 0 : for ( int i = matrix.GetHeight() - 1; i > 0; --i )
346 : : {
347 [ # # ]: 0 : for ( int j = matrix.GetWidth() - 1; j > 0; --j )
348 : : {
349 : 0 : SWNode* node = matrix(i, j);
350 : :
351 [ # # # # ]: 0 : if ( ! (node->swn_byte_assigned && ! node->swn_visited) )
352 : 0 : continue;
353 : :
354 : 0 : BroSubstring::Vec* new_al = new BroSubstring::Vec();
355 : 0 : sw_collect_single(new_al, matrix, node, params);
356 : :
357 [ # # ]: 0 : for ( vector<BroSubstring::Vec*>::iterator it = als.begin();
358 : : it != als.end(); ++it )
359 : : {
360 : 0 : BroSubstring::Vec* old_al = *it;
361 : :
362 [ # # ]: 0 : if ( old_al == 0 )
363 : 0 : continue;
364 : :
365 [ # # ]: 0 : for ( BroSubstring::VecIt it2 = old_al->begin();
366 : : it2 != old_al->end(); ++it2 )
367 : : {
368 [ # # ]: 0 : for ( BroSubstring::VecIt it3 = new_al->begin();
369 : : it3 != new_al->end(); ++it3 )
370 : : {
371 [ # # ]: 0 : if ( (*it2)->DoesCover(*it3) )
372 : : {
373 : 0 : delete_each(new_al);
374 [ # # ]: 0 : delete new_al;
375 : 0 : new_al = 0;
376 : 0 : goto end_loop;
377 : : }
378 : :
379 [ # # ]: 0 : if ( (*it3)->DoesCover(*it2) )
380 : : {
381 : 0 : delete_each(old_al);
382 [ # # ]: 0 : delete old_al;
383 : 0 : *it = 0;
384 : 0 : goto end_loop;
385 : : }
386 : : }
387 : : }
388 : : }
389 : :
390 : 0 : end_loop:
391 [ # # ]: 0 : if ( new_al )
392 : 0 : als.push_back(new_al);
393 : : }
394 : : }
395 : :
396 [ # # ]: 0 : for ( vector<BroSubstring::Vec*>::iterator it = als.begin();
397 : : it != als.end(); ++it )
398 : : {
399 : 0 : BroSubstring::Vec* al = *it;
400 : :
401 [ # # ]: 0 : if ( al == 0 )
402 : 0 : continue;
403 : :
404 [ # # ]: 0 : for ( BroSubstring::VecIt it2 = al->begin();
405 : : it2 != al->end(); ++it2 )
406 : 0 : result->push_back(*it2);
407 : :
408 [ # # ]: 0 : delete al;
409 : 0 : }
410 : 0 : }
411 : :
412 : : // The main Smith-Waterman algorithm.
413 : : //
414 : : BroSubstring::Vec* smith_waterman(const BroString* s1, const BroString* s2,
415 : 0 : SWParams& params)
416 : : {
417 : 0 : BroSubstring::Vec* result = new BroSubstring::Vec();
418 : :
419 [ # # # # ]: 0 : if ( ! s1 || s1->Len() < int(params._min_toklen) ||
[ # # ][ # # ]
[ # # ]
420 : : ! s2 || s2->Len() < int(params._min_toklen) )
421 : 0 : return result;
422 : :
423 : : // Length of both strings, plus one because SW needs
424 : : // an extra row and column.
425 : : //
426 : 0 : int i, len1 = s1->Len() + 1;
427 : 0 : int j, len2 = s2->Len() + 1;
428 : :
429 : 0 : int row = 0, col = 0;
430 : :
431 : 0 : byte_vec string1 = s1->Bytes();
432 : 0 : byte_vec string2 = s2->Bytes();
433 : :
434 : 0 : SWNodeMatrix matrix(s1, s2); // dynamic programming matrix.
435 : 0 : SWNode* node_max = 0; // pointer to the best score's node
436 : 0 : SWNode* node_br_max = 0; // pointer to lowest-right matching node
437 : :
438 : : // The highest score in the matrix, globally. We initialize to 1
439 : : // because we are only interested in real scores (initializing to
440 : : // -infty would mean 0 is larger, and would complicate the link
441 : : // structure in the matrix).
442 : : //
443 : 0 : int matrix_max = 1;
444 : 0 : int br_max_r = 0;
445 : 0 : int br_max_b = 0;
446 : :
447 : :
448 : : // Matrix initialization ----------------------------------------------
449 : :
450 : : // Assign IDs to each cell -- this is only for debugging purposes
451 : : // and can go later.
452 : :
453 : 0 : int counter = 1;
454 : :
455 [ # # ]: 0 : for ( i = 1; i < len1; ++i )
456 [ # # ]: 0 : for ( j = 1; j < len2; ++j )
457 : 0 : matrix(i, j)->id = counter++;
458 : :
459 : : // Subsequence calculation --------------------------------------------
460 : :
461 [ # # ]: 0 : for ( i = 1; i < len1; ++i )
462 : : {
463 [ # # ]: 0 : for ( j = 1; j < len2; ++j )
464 : : {
465 : : // Current node, top/left neighbours.
466 : : //
467 : 0 : SWNode* current = matrix(i, j);
468 : 0 : SWNode* node_tl = matrix(i-1, j-1);
469 : 0 : SWNode* node_l = matrix(i, j-1);
470 : 0 : SWNode* node_t = matrix(i-1, j);
471 : :
472 : : // Scores of neighbouring nodes.
473 : : //
474 : 0 : int score_t = node_t->swn_score;
475 : 0 : int score_l = node_l->swn_score;
476 : 0 : int score_tl = node_tl->swn_score;
477 : :
478 : : // If strings at current indices match, assign new
479 : : // score to current node. Minus-one adjustments
480 : : // are necessary since matrix has one extra
481 : : // row + column.
482 : : //
483 [ # # ]: 0 : if ( string1[i-1] == string2[j-1] )
484 : : {
485 : : // We have a match: improve previous score.
486 : : //
487 : 0 : score_tl += 1;
488 : :
489 : : // If we're continuing a chain of matches, rate
490 : : // higher. This favours longer consecutive
491 : : // substrings.
492 : : //
493 [ # # ]: 0 : if ( node_tl->swn_byte_assigned )
494 : 0 : score_tl += 99;
495 : :
496 : : // Store the byte we've matched in the node for
497 : : // easier access.
498 : : //
499 : 0 : current->swn_byte = string1[i-1];
500 : 0 : current->swn_byte_assigned = true;
501 : : }
502 : :
503 : : // Pick the score among the neighbours that is now highest.
504 : : // This is the core of Smith-Waterman.
505 : : //
506 [ # # ]: 0 : if ( current->swn_byte_assigned )
507 : 0 : current->swn_score = score_tl;
508 : : else
509 : 0 : current->swn_score = max(max(score_t, score_l), score_tl);
510 : :
511 : : // Establish predecessor chain according to neighbor
512 : : // with best score.
513 : : //
514 [ # # ][ # # ]: 0 : if ( current->swn_score == score_tl &&
515 : : current->swn_byte_assigned )
516 : : {
517 : : // If we had matched bytes (*and* it's the
518 : : // best neighbor), marke the node accordingly
519 : : //
520 [ # # ][ # # ]: 0 : if ( i >= br_max_b && j >= br_max_r )
521 : : {
522 : 0 : node_br_max = current;
523 : 0 : br_max_b = i;
524 : 0 : br_max_r = j;
525 : : }
526 : :
527 : 0 : current->swn_prev = node_tl;
528 : : }
529 [ # # ]: 0 : else if ( current->swn_score == score_t )
530 : 0 : current->swn_prev = node_t;
531 : : else
532 : 0 : current->swn_prev = node_l;
533 : :
534 : : // Check if we have a new global maximum -- we
535 : : // specifically track the node that is the global
536 : : // maximum so we now from where to backtrack at
537 : : // the end of the matrix iteration.
538 : : //
539 [ # # ]: 0 : if ( current->swn_score > matrix_max )
540 : : {
541 : 0 : node_max = current;
542 : 0 : matrix_max = current->swn_score;
543 : : }
544 : :
545 : : #if 0
546 : : printf("%4i/%.5i%c/%.5i[%c%c] ",
547 : : current->swn_score,
548 : : current->id,
549 : : current->swn_byte_assigned ? '*' : ' ',
550 : : current->swn_prev ? current->swn_prev->id : 0,
551 : : string1[i-1], string2[j-1]);
552 : : #endif
553 : : //printf("%.5i ", current->swn_score);
554 : : }
555 : :
556 : : #if 0
557 : : printf("\n");
558 : : #endif
559 : : }
560 : :
561 : : // Result generation.
562 : :
563 : : // How we do this depends on the mode we operate in. In SW_SINGLE, we
564 : : // follow the path from the best node until there is no predecessor
565 : : // (that is, when we hit a node in row 0), and stop. In SW_MULTIPLE,
566 : : // we collect all non-redundant common subsequences.
567 : :
568 [ # # ]: 0 : if ( params._sw_variant == SW_MULTIPLE )
569 : 0 : sw_collect_multiple(result, matrix, params);
570 : : else
571 : 0 : sw_collect_single(result, matrix, node_max, params);
572 : :
573 [ # # ]: 0 : if ( len1 > len2 )
574 : 0 : sort(result->begin(), result->end(), BroSubstringCmp(0));
575 : : else
576 : 0 : sort(result->begin(), result->end(), BroSubstringCmp(1));
577 : :
578 : 0 : return result;
579 [ + - ][ + - ]: 6 : }
|
