NVBIO
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
mapq.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 <nvBowtie/bowtie2/cuda/defs.h>
31 #include <nvBowtie/bowtie2/cuda/stats.h>
32 
33 #include <nvbio/io/alignments.h>
34 #include <nvbio/io/output/output_types.h>
35 
36 namespace nvbio {
37 namespace bowtie2 {
38 namespace cuda {
39 
42 template <typename ScoringSchemeType>
43 struct BowtieMapq3
44 {
45  typedef ScoringSchemeType scoring_scheme_type;
46 
49  BowtieMapq3(const ScoringSchemeType sc) : m_sc(sc) {}
50 
53  NVBIO_HOST_DEVICE
54  uint32 operator() (
55  const io::BestPairedAlignments& best_alignments,
56  const uint32 read_len,
57  const uint32 o_read_len) const
58  {
59  // no second best alignment, perfect score
60  const int unpaired_one_perfect = 44;
61 
62  // no second best alignment
63  const int unpaired_one[11] =
64  {
65  43, 42, 41, 36, 32, 27, 20, 11, 4, 1, 0
66  };
67 
68  // two alignment scores, the best of which has perfect score
69  const int unpaired_two_perfect[11] =
70  {
71  2, 16, 23, 30, 31, 32, 34, 36, 38, 40, 42
72  };
73 
74  // two alignment scores, the best of which has a non perfect score
75  const int unpaired_two[11][11] =
76  {
77  { 2, 2, 2, 1, 1, 0, 0, 0, 0, 0, 0 },
78  { 20, 14, 7, 3, 2, 1, 0, 0, 0, 0, 0 },
79  { 20, 16, 10, 6, 3, 1, 0, 0, 0, 0, 0 },
80  { 20, 17, 13, 9, 3, 1, 1, 0, 0, 0, 0 },
81  { 21, 19, 15, 9, 5, 2, 2, 0, 0, 0, 0 },
82  { 22, 21, 16, 11, 10, 5, 0, 0, 0, 0, 0 },
83  { 23, 22, 19, 16, 11, 0, 0, 0, 0, 0, 0 },
84  { 24, 25, 21, 30, 0, 0, 0, 0, 0, 0, 0 },
85  { 30, 26, 29, 0, 0, 0, 0, 0, 0, 0, 0 },
86  { 30, 27, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
87  { 30, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
88  };
89 
90  // paired alignment, one perfect score
91  const int paired_one_perfect = 44;
92 
93  //
94  // TODO: if paired, return paired_one_perfect
95  //
96  if (best_alignments.is_paired())
97  return paired_one_perfect;
98 
99  const float max_score = (float)m_sc.perfect_score( read_len );
100  const float min_score = (float)m_sc.min_score( read_len );
101  const float norm_factor = 10.0f / (max_score - min_score);
102 
103  // check whether the best score is beyond the minimum threshold
104  if (best_alignments.best_score() < min_score)
105  return 0;
106 
107  const int best = std::max( (int)max_score - best_alignments.best_score(), 0 ); // negated best score
108  const int best_bin = (int)(float(best) * norm_factor + 0.5f);
109 
110  const bool has_second = best_alignments.has_second();
111 
112 // if (info)
113 // {
114 // info->best = best;
115 // info->second_best = has_second ? best_alignments.second_score() : 255;
116 // }
117 
118  // check whether there are two alignment scores or just one
119  if (has_second)
120  {
121  const int diff = best_alignments.best_score() - best_alignments.second_score();
122  const int diff_bin = (int)(float(diff) * norm_factor + 0.5f);
123 
124  if (best == max_score)
125  return unpaired_two_perfect[best_bin];
126  else
127  return unpaired_two[diff_bin][best_bin];
128  }
129  else
130  {
131  if (best == max_score)
132  return unpaired_one_perfect;
133  else
134  return unpaired_one[best_bin];
135  }
136  }
137 
138 private:
139  ScoringSchemeType m_sc;
140 };
141 
144 template <typename ScoringSchemeType>
145 struct BowtieMapq2
146 {
147  typedef ScoringSchemeType scoring_scheme_type;
148 
151  BowtieMapq2(const ScoringSchemeType sc) : m_sc(sc) {}
152 
155  NVBIO_HOST_DEVICE
156  uint32 operator() (
157  const io::BestPairedAlignments& best_alignments,
158  const uint32 read_len,
159  const uint32 o_read_len) const
160  {
161  bool has_second = best_alignments.has_second();
162 
163  const float max_score = (float)m_sc.perfect_score( read_len ) +
164  (best_alignments.is_paired() ? m_sc.perfect_score( o_read_len ) : 0);
165  const float min_score = (float)m_sc.min_score( read_len ) +
166  (best_alignments.is_paired() ? m_sc.min_score( o_read_len ) : 0);
167 
168  const float diff = (max_score - min_score); // range of scores
169 
170  const float best = (float)best_alignments.best_score();
171  float second_best = min_score-1;
172 
173  if (best < min_score)
174  return 0;
175 
176  // best score but normalized so that 0 = worst valid score
177  const float best_over = best - min_score;
178 
179 // if (info)
180 // {
181 // info->best = (int32)best;
182 // info->second_best = (int32)second_best;
183 // }
184 
185  if (m_sc.m_monotone)
186  {
187  // global alignment
188  if (!has_second)
189  {
190  if (best_over >= diff * 0.8f) return 42;
191  else if (best_over >= diff * 0.7f) return 40;
192  else if (best_over >= diff * 0.6f) return 24;
193  else if (best_over >= diff * 0.5f) return 23;
194  else if (best_over >= diff * 0.4f) return 8;
195  else if (best_over >= diff * 0.3f) return 3;
196  else return 0;
197  }
198  else
199  {
200  second_best = (float)best_alignments.second_score();
201 // if (info)
202 // info->second_best = (int32)second_best;
203 
204  const float best_diff = fabsf(fabsf( best ) - fabsf( second_best ));
205 
206  if (best_diff >= diff * 0.9f)
207  return (best_over == diff) ? 39 : 33;
208  else if (best_diff >= diff * 0.8f)
209  return (best_over == diff) ? 38 : 27;
210  else if (best_diff >= diff * 0.7f)
211  return (best_over == diff) ? 37 : 26;
212  else if (best_diff >= diff * 0.6f)
213  return (best_over == diff) ? 36 : 22;
214  else if (best_diff >= diff * 0.5f)
215  {
216  // top third is still pretty good
217  if (best_over == diff) return 35;
218  else if (best_over >= diff * 0.84f) return 25;
219  else if (best_over >= diff * 0.68f) return 16;
220  else return 5;
221  }
222  else if (best_diff >= diff * 0.4f)
223  {
224  // top third is still pretty good
225  if (best_over == diff) return 34;
226  else if (best_over >= diff * 0.84f) return 21;
227  else if (best_over >= diff * 0.68f) return 14;
228  else return 4;
229  }
230  else if (best_diff >= diff * 0.3f)
231  {
232  // top third is still pretty good
233  if (best_over == diff) return 32;
234  else if (best_over >= diff * 0.88f) return 18;
235  else if (best_over >= diff * 0.67f) return 15;
236  else return 3;
237  }
238  else if (best_diff >= diff * 0.2f)
239  {
240  // top third is still pretty good
241  if (best_over == diff) return 31;
242  else if (best_over >= diff * 0.88f) return 17;
243  else if (best_over >= diff * 0.67f) return 11;
244  else return 0;
245  }
246  else if (best_diff >= diff * 0.1f)
247  {
248  // top third is still pretty good
249  if (best_over == diff) return 30;
250  else if (best_over >= diff * 0.88f) return 12;
251  else if (best_over >= diff * 0.67f) return 7;
252  else return 0;
253  }
254  else if (best_diff > 0)
255  {
256  // top third is still pretty good
257  return (best_over >= diff * 0.67f) ? 6 : 2;
258  }
259  else
260  {
261  // top third is still pretty good
262  return (best_over >= diff * 0.67f) ? 1 : 0;
263  }
264  }
265  }
266  else
267  {
268  // local alignment
269  if (!has_second)
270  {
271  if (best_over >= diff * 0.8f) return 44;
272  else if (best_over >= diff * 0.7f) return 42;
273  else if (best_over >= diff * 0.6f) return 41;
274  else if (best_over >= diff * 0.5f) return 36;
275  else if (best_over >= diff * 0.4f) return 28;
276  else if (best_over >= diff * 0.3f) return 24;
277  else return 22;
278  }
279  else
280  {
281  second_best = (float)best_alignments.second_score();
282 // if (info)
283 // info->second_best = (int32)second_best;
284 
285  const float best_diff = fabsf(fabsf( best ) - fabsf( second_best ));
286 
287  if (best_diff >= diff * 0.9f) return 40;
288  else if (best_diff >= diff * 0.8f) return 39;
289  else if (best_diff >= diff * 0.7f) return 38;
290  else if (best_diff >= diff * 0.6f) return 37;
291  else if (best_diff >= diff * 0.5f)
292  {
293  if (best_over == diff) return 35;
294  else if (best_over >= diff * 0.50f) return 25;
295  else return 20;
296  }
297  else if (best_diff >= diff * 0.4f)
298  {
299  if (best_over == diff) return 34;
300  else if (best_over >= diff * 0.50f) return 21;
301  else return 19;
302  }
303  else if (best_diff >= diff * 0.3f)
304  {
305  if (best_over == diff) return 33;
306  else if (best_over >= diff * 0.5f) return 18;
307  else return 16;
308  }
309  else if (best_diff >= diff * 0.2f)
310  {
311  if (best_over == diff) return 32;
312  else if (best_over >= diff * 0.5f) return 17;
313  else return 12;
314  }
315  else if (best_diff >= diff * 0.1f)
316  {
317  if (best_over == diff) return 31;
318  else if (best_over >= diff * 0.5f) return 14;
319  else return 9;
320  }
321  else if (best_diff > 0)
322  return (best_over >= diff * 0.5f) ? 11 : 2;
323  else
324  return (best_over >= diff * 0.5f) ? 1 : 0;
325  }
326  }
327  }
328 
329 private:
330  ScoringSchemeType m_sc;
331 };
332 
333 } // namespace cuda
334 } // namespace bowtie2
335 } // namespace nvbio
Generated on Wed Feb 25 2015 08:33:05 for NVBIO by   doxygen 1.8.4