NVBIO
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
sw_banded_inl.h
Go to the documentation of this file.
1 /*
2  * nvbio
3  * Copyright (c) 2011-2014, NVIDIA CORPORATION. All rights reserved.
4  *
5  * Redistribution and use in source and binary forms, with or without
6  * modification, are permitted provided that the following conditions are met:
7  * * Redistributions of source code must retain the above copyright
8  * notice, this list of conditions and the following disclaimer.
9  * * Redistributions in binary form must reproduce the above copyright
10  * notice, this list of conditions and the following disclaimer in the
11  * documentation and/or other materials provided with the distribution.
12  * * Neither the name of the NVIDIA CORPORATION nor the
13  * names of its contributors may be used to endorse or promote products
14  * derived from this software without specific prior written permission.
15  *
16  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
17  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
18  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
19  * DISCLAIMED. IN NO EVENT SHALL NVIDIA CORPORATION BE LIABLE FOR ANY
20  * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
21  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
22  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
23  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
25  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26  */
27 
28 #pragma once
29 
30 #include <nvbio/basic/types.h>
31 #include <nvbio/basic/numbers.h>
32 
33 namespace nvbio {
34 namespace aln {
35 
36 namespace priv {
37 namespace banded {
38 
41 
42 // initialize the zero-th row of the banded DP-matrix with Smith-Waterman scoring
43 //
44 template <uint32 BAND_LEN, AlignmentType TYPE, typename band_type, typename scoring_type, typename score_type>
45 NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
46 void init_row_zero(
47  band_type& band,
48  const scoring_type& scoring,
49  const score_type infimum)
50 {
51  #pragma unroll
52  for (uint32 j = 0; j < BAND_LEN; ++j)
53  band[j] = TYPE == GLOBAL ? j * scoring.deletion() : 0;
54 }
55 
64 template <uint32 BAND_LEN, AlignmentType TYPE>
65 struct SmithWatermanScoringContext
66 {
67  template <typename band_type, typename scoring_type, typename score_type>
68  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
69  void init(
70  const uint32 i,
71  band_type& band,
72  const scoring_type& scoring,
73  const score_type infimum)
74  {
75  init_row_zero<BAND_LEN,TYPE>( band, scoring, infimum );
76  }
77 
78  template <typename band_type>
79  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
80  void previous_row(
81  const uint32 i,
82  const band_type& band) {}
83 
84  template <typename band_type>
85  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
86  void last_row(
87  const uint32 i,
88  const band_type& band) {}
89 
90  template <typename score_type>
91  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
92  void new_cell(
93  const uint32 i,
94  const uint32 j,
95  const score_type score,
96  const DirectionVector dir) {}
97 };
98 
103 template <uint32 BAND_LEN, AlignmentType TYPE, typename checkpoint_type>
104 struct SmithWatermanCheckpointedScoringContext
105 {
106  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
107  SmithWatermanCheckpointedScoringContext(checkpoint_type checkpoints) :
108  m_checkpoints( checkpoints ) {}
109 
110  template <typename band_type, typename scoring_type, typename score_type>
111  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
112  void init(
113  const uint32 i,
114  band_type& band,
115  const scoring_type& scoring,
116  const score_type infimum)
117  {
118  // check whether this is the first row
119  if (i == 0)
120  init_row_zero<BAND_LEN,TYPE>( band, scoring, infimum );
121  else
122  {
123  // load from checkpoint
124  #pragma unroll
125  for (uint32 j = 0; j < BAND_LEN; ++j)
126  band[j] = m_checkpoints[j];
127  }
128  }
129 
130  template <typename band_type>
131  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
132  void previous_row(
133  const uint32 i,
134  const band_type& band) {}
135 
136  template <typename band_type>
137  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
138  void last_row(
139  const uint32 i,
140  const band_type& band)
141  {
142  const short infimum = Field_traits<short>::min() + 32;
143 
144  // save the last row
145  #pragma unroll
146  for (uint32 j = 0; j < BAND_LEN; ++j)
147  m_checkpoints[j] = nvbio::max( band[j], infimum );
148  }
149 
150  // dir is the backtracking arrow for the H matrix at this cell
151  template <typename score_type>
152  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
153  void new_cell(
154  const uint32 i,
155  const uint32 j,
156  const score_type score,
157  const DirectionVector dir) {}
158 
159  checkpoint_type m_checkpoints;
160 };
161 
166 template <uint32 BAND_LEN, AlignmentType TYPE, uint32 CHECKPOINTS, typename checkpoint_type>
167 struct SmithWatermanCheckpointContext
168 {
169  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
170  SmithWatermanCheckpointContext(checkpoint_type checkpoints) :
171  m_checkpoints( checkpoints ) {}
172 
173  template <typename band_type, typename scoring_type, typename score_type>
174  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
175  void init(
176  const uint32 i,
177  band_type& band,
178  const scoring_type& scoring,
179  const score_type infimum)
180  {
181  init_row_zero<BAND_LEN,TYPE>( band, scoring, infimum );
182  }
183 
184  template <typename band_type>
185  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
186  void previous_row(
187  const uint32 i,
188  const band_type& band)
189  {
190  // save checkpoint
191  if ((i & (CHECKPOINTS-1)) == 0u)
192  {
193  const short infimum = Field_traits<short>::min() + 32;
194 
195  for (uint32 j = 0; j < BAND_LEN; ++j)
196  m_checkpoints[BAND_LEN*(i/CHECKPOINTS) + j] = nvbio::max( band[j], infimum );
197  }
198  }
199 
200  template <typename band_type>
201  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
202  void last_row(
203  const uint32 i,
204  const band_type& band)
205  {
206  const short infimum = Field_traits<short>::min() + 32;
207 
208  // save the last row
209  const uint32 checkpoint_row = (i+CHECKPOINTS-1)/CHECKPOINTS;
210  for (uint32 j = 0; j < BAND_LEN; ++j)
211  m_checkpoints[BAND_LEN*checkpoint_row + j] = nvbio::max( band[j], infimum );
212  }
213 
214  // dir is the backtracking arrow for the H matrix at this cell
215  template <typename score_type>
216  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
217  void new_cell(
218  const uint32 i,
219  const uint32 j,
220  const score_type score,
221  const DirectionVector dir) {}
222 
223  checkpoint_type m_checkpoints;
224 };
225 
230 template <uint32 BAND_LEN, AlignmentType TYPE, uint32 CHECKPOINTS, typename checkpoint_type, typename submatrix_type>
231 struct SmithWatermanSubmatrixContext
232 {
233  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
234  SmithWatermanSubmatrixContext(
235  const checkpoint_type checkpoints,
236  const uint32 checkpoint_id,
237  submatrix_type submatrix) :
238  m_checkpoints( checkpoints ),
239  m_checkpoint_id( checkpoint_id ),
240  m_submatrix( submatrix ) {}
241 
242  template <typename band_type, typename scoring_type, typename score_type>
243  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
244  void init(
245  const uint32 i,
246  band_type& band,
247  const scoring_type& scoring,
248  const score_type infimum)
249  {
250  // restore the checkpoint
251  #pragma unroll
252  for (uint32 j = 0; j < BAND_LEN; ++j)
253  band[j] = m_checkpoints[m_checkpoint_id*BAND_LEN + j];
254  }
255 
256  template <typename band_type>
257  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
258  void previous_row(
259  const uint32 i,
260  const band_type& band) {}
261 
262  template <typename band_type>
263  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
264  void last_row(
265  const uint32 i,
266  const band_type& band) {}
267 
268  template <typename score_type>
269  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
270  void new_cell(
271  const uint32 i,
272  const uint32 j,
273  const score_type score,
274  const DirectionVector dir)
275  {
276  // save the direction vectors for H,E
277  const uint32 offset = m_checkpoint_id * CHECKPOINTS;
278  m_submatrix[ (i - offset)*BAND_LEN + j ] = dir;
279  }
280 
281  checkpoint_type m_checkpoints;
282  uint32 m_checkpoint_id;
283  submatrix_type m_submatrix;
284 };
285 
289 template <uint32 BAND_LEN, AlignmentType TYPE>
290 struct sw_alignment_score_dispatch
291 {
340  template <
341  typename pattern_type,
342  typename qual_type,
343  typename text_type,
344  typename scoring_type,
345  typename context_type,
346  typename sink_type>
347  NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
348  static
349  bool run(
350  const scoring_type& scoring,
351  pattern_type pattern,
352  qual_type quals,
353  text_type text,
354  const uint32 window_begin,
355  const uint32 window_end,
356  const uint32 pos,
357  const int32 min_score,
358  context_type& context,
359  sink_type& sink)
360  {
361  const uint32 pattern_len = pattern.length();
362  const uint32 text_len = text.length();
363  const uint32 start = pos;
364 
365  if (text_len < pattern_len)
366  return false;
367 
368  typedef int32 score_type;
369 
370  typedef typename Reference_cache<BAND_LEN>::type text_cache_type;
371  uint32 text_cache_storage[Reference_cache<BAND_LEN>::BAND_WORDS];
372  text_cache_type text_cache( &text_cache_storage[0] );
373 
374  // load first band of text
375  for (uint32 j = 0; j < BAND_LEN-1; ++j)
376  text_cache[j] = text[start+window_begin+j];
377 
378  const score_type G = scoring.deletion();
379  const score_type I = scoring.insertion();
380  const score_type infimum = Field_traits<short>::min() + nvbio::max( G, I );
381 
382  // initialize the first band (corresponding to the 0-th row of the DP matrix)
383  score_type band[BAND_LEN];
384 
385  // initialize band
386  context.init(
387  window_begin,
388  band,
389  scoring,
390  infimum );
391 
392  // loop across the short edge of the DP matrix: each band is a segment of the long columns
393  for (uint32 i = window_begin; i < window_end; ++i)
394  {
395  // pass the previous row to the context
396  context.previous_row( i, band );
397 
398  // load the new pattern character
399  const uint8 q = pattern[i];
400  const uint8 qq = quals[i];
401 
402  const score_type V = scoring.match(qq);
403 
404  score_type top, left, diagonal, hi;
405 
406  // j == 0 case
407  {
408  const uint8 g = text_cache[0];
409  const score_type S_ij = (g == q) ? V : scoring.mismatch( g, q, qq );
410  diagonal = band[0] + S_ij;
411  top = band[1] + G;
412  hi = nvbio::max( top, diagonal );
413  if (TYPE == LOCAL)
414  {
415  hi = nvbio::max( hi, score_type(0) ); // clamp to zero
416  sink.report( hi, make_uint2( i+1, i+1 ) );
417  }
418  band[0] = hi;
419 
420  // pass the new cell to the context
421  context.new_cell(
422  i, 0u,
423  hi,
424  top > diagonal ? INSERTION : SUBSTITUTION );
425  }
426 
427  #pragma unroll
428  for (uint32 j = 1; j < BAND_LEN-1; ++j)
429  {
430  const uint8 g = text_cache[j]; text_cache[j-1] = g;
431  const score_type S_ij = (g == q) ? V : scoring.mismatch( g, q, qq );
432  diagonal = band[j] + S_ij;
433  top = band[j+1] + G;
434  left = band[j-1] + I;
435  hi = nvbio::max3( top, left, diagonal );
436  if (TYPE == LOCAL)
437  {
438  hi = nvbio::max( hi, score_type(0) ); // clamp to zero
439  sink.report( hi, make_uint2( i + j + 1, i+1 ) );
440  }
441  band[j] = hi;
442 
443  // pass the new cell to the context
444  context.new_cell(
445  i, j,
446  hi,
447  top > left ?
448  (top > diagonal ? INSERTION : SUBSTITUTION) :
449  (left > diagonal ? DELETION : SUBSTITUTION) );
450  }
451 
452  // load the new text character and fill last entry of the band cache
453  const uint8 g = start+i+BAND_LEN-1 < text_len ? text[start+i+BAND_LEN-1] : 255u;
454  text_cache[ BAND_LEN-2 ] = g;
455 
456  // j == BAND_LEN-1 case
457  {
458  const score_type S_ij = (g == q) ? V : scoring.mismatch( g, q, qq );
459  diagonal = band[BAND_LEN-1] + S_ij;
460  left = band[BAND_LEN-2] + I;
461  hi = nvbio::max( left, diagonal );
462  if (TYPE == LOCAL)
463  {
464  hi = nvbio::max( hi, score_type(0) ); // clamp to zero
465  sink.report( hi, make_uint2( i + BAND_LEN, i+1 ) );
466  }
467  band[ BAND_LEN-1 ] = hi;
468 
469  // pass the new cell to the context
470  context.new_cell(
471  i, BAND_LEN-1,
472  hi,
473  left > diagonal ? DELETION : SUBSTITUTION );
474  }
475  }
476 
477  if (window_end < pattern_len)
478  {
479  // compute the maximum score we got
480  score_type max_score = band[0];
481  #pragma unroll
482  for (uint32 j = 1; j < BAND_LEN; ++j)
483  max_score = nvbio::max( max_score, band[j] );
484 
485  // and check whether min_score is within reach
486  const score_type threshold_score = score_type( min_score ) + score_type(pattern_len - window_end)*scoring.match(0);
487  if (max_score < threshold_score)
488  return false;
489  }
490 
491  // pass the last row to the context
492  context.last_row( window_end, band );
493 
494  if (window_end == pattern_len)
495  {
496  if (TYPE == GLOBAL)
497  sink.report( band[BAND_LEN-1], make_uint2( pattern_len + BAND_LEN-1, pattern_len ) );
498  else if (TYPE == SEMI_GLOBAL)
499  {
500  const uint32 m = nvbio::min( pattern_len + BAND_LEN - 1u, text_len ) - (pattern_len-1u);
501 
502  // get the highest score along the long edge of the path graph
503  sink.report( band[0], make_uint2( pattern_len + 0, pattern_len ) );
504  #pragma unroll
505  for (uint32 j = 1; j < BAND_LEN; ++j)
506  {
507  if (j < m)
508  sink.report( band[j], make_uint2( pattern_len + j, pattern_len ) );
509  }
510  }
511  }
512  return true;
513  }
514 };
515 
517 
518 } // namespace banded
519 
522 
535 template <
536  uint32 BAND_LEN,
537  AlignmentType TYPE,
538  typename scoring_type,
539  typename pattern_type,
540  typename qual_type,
541  typename text_type,
542  typename sink_type>
543 NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
544 bool banded_alignment_score(
545  const SmithWatermanAligner<TYPE,scoring_type>& aligner,
546  pattern_type pattern,
547  qual_type quals,
548  text_type text,
549  const int32 min_score,
550  sink_type& sink)
551 {
552  priv::banded::SmithWatermanScoringContext<BAND_LEN,TYPE> context;
553 
554  return priv::banded::sw_alignment_score_dispatch<BAND_LEN,TYPE>::run( aligner.scheme, pattern, quals, text, 0u, pattern.length(), 0u, min_score, context, sink );
555 }
556 
570 template <
571  uint32 BAND_LEN,
572  AlignmentType TYPE,
573  typename scoring_type,
574  typename pattern_type,
575  typename qual_type,
576  typename text_type,
577  typename sink_type,
578  typename checkpoint_type>
579 NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
580 bool banded_alignment_score(
581  const SmithWatermanAligner<TYPE,scoring_type>& aligner,
582  pattern_type pattern,
583  qual_type quals,
584  text_type text,
585  const int32 min_score,
586  const uint32 window_begin,
587  const uint32 window_end,
588  sink_type& sink,
589  checkpoint_type checkpoint)
590 {
591  priv::banded::SmithWatermanCheckpointedScoringContext<BAND_LEN,TYPE,checkpoint_type> context( checkpoint );
592 
593  return priv::banded::sw_alignment_score_dispatch<BAND_LEN,TYPE>::run( aligner.scheme, pattern, quals, text, window_begin, window_end, 0u, min_score, context, sink );
594 }
595 
618 template <
619  uint32 BAND_LEN,
620  uint32 CHECKPOINTS,
621  AlignmentType TYPE,
622  typename scoring_type,
623  typename pattern_type,
624  typename qual_type,
625  typename text_type,
626  typename sink_type,
627  typename checkpoint_type>
628 NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
629 bool banded_alignment_checkpoints(
630  const SmithWatermanAligner<TYPE,scoring_type>& aligner,
631  pattern_type pattern,
632  qual_type quals,
633  text_type text,
634  const int32 min_score,
635  sink_type& sink,
636  checkpoint_type checkpoints)
637 {
638  priv::banded::SmithWatermanCheckpointContext<BAND_LEN,TYPE,CHECKPOINTS,checkpoint_type> context( checkpoints );
639 
640  return priv::banded::sw_alignment_score_dispatch<BAND_LEN,TYPE>::run( aligner.scheme, pattern, quals, text, 0u, pattern.length(), 0u, min_score, context, sink );
641 }
642 
675 template <
676  uint32 BAND_LEN,
677  uint32 CHECKPOINTS,
678  AlignmentType TYPE,
679  typename scoring_type,
680  typename pattern_string,
681  typename qual_string,
682  typename text_string,
683  typename checkpoint_type,
684  typename submatrix_type>
685 NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
686 uint32 banded_alignment_submatrix(
687  const SmithWatermanAligner<TYPE,scoring_type>& aligner,
688  pattern_string pattern,
689  qual_string quals,
690  text_string text,
691  const int32 min_score,
692  checkpoint_type checkpoints,
693  const uint32 checkpoint_id,
694  submatrix_type submatrix)
695 {
696  priv::banded::SmithWatermanSubmatrixContext<BAND_LEN,TYPE,CHECKPOINTS,checkpoint_type,submatrix_type> context( checkpoints, checkpoint_id, submatrix );
697  const uint32 window_begin = checkpoint_id * CHECKPOINTS;
698  const uint32 window_end = nvbio::min( window_begin + CHECKPOINTS, uint32(pattern.length()) );
699  NullSink sink;
700 
701  priv::banded::sw_alignment_score_dispatch<BAND_LEN,TYPE>::run( aligner.scheme, pattern, quals, text, window_begin, window_end, 0u, min_score, context, sink );
702 
703  return window_end - window_begin;
704 }
705 
740 template <
741  uint32 BAND_LEN,
742  uint32 CHECKPOINTS,
743  AlignmentType TYPE,
744  typename scoring_type,
745  typename checkpoint_type,
746  typename submatrix_type,
747  typename backtracer_type>
748 NVBIO_FORCEINLINE NVBIO_HOST_DEVICE
749 bool banded_alignment_traceback(
750  const SmithWatermanAligner<TYPE,scoring_type>& aligner,
751  checkpoint_type checkpoints,
752  const uint32 checkpoint_id,
753  submatrix_type submatrix,
754  const uint32 submatrix_height,
755  uint8& state,
756  uint2& sink,
757  backtracer_type& backtracer)
758 {
759  // Backtrack from the second checkpoint to the first looking up the flow submatrices
760 
761  int32 current_entry = sink.x - sink.y;
762  int32 current_row = sink.y - checkpoint_id*CHECKPOINTS - 1u;
763 
764  NVBIO_CUDA_DEBUG_ASSERT( current_entry >= 0 &&
765  current_entry < BAND_LEN, "sw::banded_alignment_backtrack(): sink (%u,%u) -> local x coordinate %d not in [0,%d[\n", sink.x, sink.y, current_entry, BAND_LEN );
766  NVBIO_CUDA_DEBUG_ASSERT( current_row >= 0, "sw::banded_alignment_backtrack(): sink (%u,%u) -> local y coordinate %d not in [0,%u[\n", sink.x, sink.y, current_row, submatrix_height );
767  NVBIO_CUDA_DEBUG_ASSERT( current_row < submatrix_height, "sw::banded_alignment_backtrack(): sink (%u,%u) -> local y coordinate %d not in [0,%u[\n", sink.x, sink.y, current_row, submatrix_height );
768 
769  while (current_row >= 0)
770  {
771  const uint32 submatrix_cell = current_row * BAND_LEN + current_entry;
772  const uint8 op = submatrix[ submatrix_cell ];
773 
774  if (TYPE == LOCAL)
775  {
776  if (op == SINK)
777  {
778  sink.y = current_row + checkpoint_id*CHECKPOINTS + 1u;
779  sink.x = current_entry + sink.y;
780  return true;
781  }
782  }
783 
784  if (op == DELETION)
785  {
786  --current_entry;
787  backtracer.push( DELETION );
788  }
789  else if (op == INSERTION)
790  {
791  ++current_entry; --current_row;
792  backtracer.push( INSERTION );
793  }
794  else
795  {
796  --current_row;
797  backtracer.push( SUBSTITUTION );
798  }
799 
800  NVBIO_CUDA_ASSERT( current_entry >= 0 && current_entry < BAND_LEN );
801  }
802  sink.y = checkpoint_id*CHECKPOINTS;
803  sink.x = current_entry + sink.y;
804  return false;
805 }
806 
808 
809 } // namespace priv
810 
811 } // namespace aln
812 } // namespace nvbio
813 
Generated on Wed Feb 25 2015 08:32:53 for NVBIO by   doxygen 1.8.4