NVBIO
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
score_paired_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 
32 
33 namespace nvbio {
34 namespace bowtie2 {
35 namespace cuda {
36 
37 namespace detail {
38 
41 
44 
47 
53 template <typename AlignerType, typename PipelineType>
54 struct BestAnchorScoreStream : public AlignmentStreamBase<SCORE_STREAM,AlignerType,PipelineType>
55 {
57  typedef typename base_type::context_type context_type;
59 
71  const uint32 _band_len,
72  const PipelineType _pipeline,
73  const AlignerType _aligner,
74  const ParamsPOD _params) :
75  base_type( _pipeline, _aligner, _params ), m_band_len( _band_len ) {}
76 
80  uint32 max_pattern_length() const { return base_type::m_pipeline.reads.max_read_len(); }
81 
85  uint32 max_text_length() const { return base_type::m_pipeline.reads.max_read_len() + m_band_len; }
86 
90  uint32 size() const { return base_type::m_pipeline.hits_queue_size; }
91 
96  const uint32 i,
97  context_type* context) const
98  {
99  context->idx = base_type::m_pipeline.idx_queue[i];
100 
101  // initialize the sink
102  context->sink = aln::BestSink<int32>();
103 
104  // fetch the hit to process
105  HitReference<HitQueuesDeviceView> hit = base_type::m_pipeline.scoring_queues.hits[ context->idx ];
106 
107  const uint32 read_rc = hit.seed.rc;
108  const uint32 read_id = hit.read_id;
109 
111  io::BestAlignments( base_type::m_pipeline.best_alignments[ read_id ], base_type::m_pipeline.best_alignments[ read_id + base_type::m_pipeline.best_stride ] ),
112  io::BestAlignments( base_type::m_pipeline.best_alignments_o[ read_id ], base_type::m_pipeline.best_alignments_o[ read_id + base_type::m_pipeline.best_stride ] ) );
113 
114  // compute the optimal score of the opposite mate
115  const uint2 o_read_range = base_type::m_pipeline.reads_o.get_range( read_id );
116  const uint32 o_read_len = o_read_range.y - o_read_range.x;
117  const int32 o_optimal_score = base_type::m_pipeline.scoring_scheme.perfect_score( o_read_len );
118  const int32 o_worst_score = base_type::m_pipeline.scoring_scheme.min_score( o_read_len );
119 
120  const uint2 a_read_range = base_type::m_pipeline.reads.get_range( read_id );
121  const uint32 a_read_len = a_read_range.y - a_read_range.x;
122  const int32 a_optimal_score = base_type::m_pipeline.scoring_scheme.perfect_score( a_read_len );
123  const int32 a_worst_score = base_type::m_pipeline.scoring_scheme.min_score( a_read_len );
124 
125  const int32 target_pair_score = nvbio::min(
126  compute_target_score( best, a_worst_score, o_worst_score ) + 1,
127  a_optimal_score + o_optimal_score );
128 
129  int32 target_mate_score = target_pair_score - o_optimal_score;
130 
131  #if 1
132  //
133  // bound the score of this mate by the worst score allowed for its own read length.
134  // NOTE: disabling this is equivalent to allowing a worst score proportional to the total
135  // length of the read pair.
136  //
137 
138  target_mate_score = nvbio::max( target_mate_score, a_worst_score );
139  #endif
140 
141  // setup the read info
142  context->mate = 0;
143  context->read_rc = read_rc;
144  context->read_id = read_id;
145  context->read_range = a_read_range;
146 
147  // setup the genome range
148  const uint32 g_pos = hit.loc;
149 
150  context->genome_begin = g_pos > m_band_len/2 ? g_pos - m_band_len/2 : 0u;
151  context->genome_end = nvbio::min( context->genome_begin + m_band_len + a_read_len, base_type::m_pipeline.genome_length );
152 
153  // skip locations that we have already visited
154  const uint32 mate = base_type::m_pipeline.anchor;
155  const bool skip = (mate == best.m_a1.m_mate && (read_rc == best.m_a1.m_rc && g_pos == best.m_a1.m_align)) ||
156  (mate == best.m_o1.m_mate && (read_rc == best.m_o1.m_rc && g_pos == best.m_o1.m_align)) ||
157  (mate == best.m_a2.m_mate && (read_rc == best.m_a2.m_rc && g_pos == best.m_a2.m_align)) ||
158  (mate == best.m_o2.m_mate && (read_rc == best.m_o2.m_rc && g_pos == best.m_o2.m_align)) ||
159  (context->min_score > a_optimal_score);
160 
161  // setup the minimum score
162  context->min_score = skip ? Field_traits<int32>::max() : nvbio::max( target_mate_score, base_type::m_pipeline.score_limit );
163 
164  return !skip;
165  }
166 
170  void output(
171  const uint32 i,
172  const context_type* context) const
173  {
174  // write the final hit.score and hit.sink attributes
175  HitReference<HitQueuesDeviceView> hit = base_type::m_pipeline.scoring_queues.hits[ context->idx ];
176 
177  const aln::BestSink<int32> sink = context->sink;
178  hit.score = (sink.score >= context->min_score) ? sink.score : scheme_type::worst_score;
179  hit.sink = context->genome_begin + sink.sink.x;
180  // TODO: rewrite hit.loc
181  // hit.loc = context->genome_begin;
182  NVBIO_CUDA_DEBUG_PRINT_IF( base_type::m_params.debug.show_score( context->read_id, (sink.score >= context->min_score) ), "score anchor: %d (min[%d], mate[%u], rc[%u], pos[%u], [qid %u])\n", sink.score, context->min_score, base_type::m_pipeline.anchor, context->read_rc, context->genome_begin, i );
183  }
184 
186 };
187 
191 template <typename aligner_type, typename pipeline_type>
193  const uint32 band_len,
194  const pipeline_type& pipeline,
195  const aligner_type aligner,
196  const ParamsPOD params)
197 {
198  const uint32 static_band_len =
199  (band_len < 4) ? 3u :
200  (band_len < 8) ? 7u :
201  (band_len < 16) ? 15u :
202  31u;
203 
205 
206  stream_type stream(
207  static_band_len,
208  pipeline,
209  aligner,
210  params );
211 
212  typedef aln::DeviceThreadScheduler scheduler_type;
213  //typedef aln::DeviceStagedThreadScheduler scheduler_type;
214 
215  if (band_len < 4)
216  {
218 
219  batch.enact( stream, pipeline.dp_buffer_size, pipeline.dp_buffer );
220  }
221  else if (band_len < 8)
222  {
224 
225  batch.enact( stream, pipeline.dp_buffer_size, pipeline.dp_buffer );
226  }
227  else if (band_len < 16)
228  {
230 
231  batch.enact( stream, pipeline.dp_buffer_size, pipeline.dp_buffer );
232  }
233  else
234  {
236 
237  batch.enact( stream, pipeline.dp_buffer_size, pipeline.dp_buffer );
238  }
239 }
240 
244 
245 } // namespace detail
246 
247 //
248 // execute a batch of banded-alignment score calculations for the anchor mates, best mapping
249 //
250 template <typename scheme_type>
252  const uint32 band_len,
254  const ParamsPOD params)
255 {
256  if (params.alignment_type == LocalAlignment)
257  {
259  band_len,
260  pipeline,
261  pipeline.scoring_scheme.local_aligner(),
262  params );
263  }
264  else
265  {
267  band_len,
268  pipeline,
269  pipeline.scoring_scheme.end_to_end_aligner(),
270  params );
271  }
272 }
273 
274 } // namespace cuda
275 } // namespace bowtie2
276 } // namespace nvbio