OpenMS  2.4.0
SpectrumAlignment.h
Go to the documentation of this file.
1 // --------------------------------------------------------------------------
2 // OpenMS -- Open-Source Mass Spectrometry
3 // --------------------------------------------------------------------------
4 // Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
5 // ETH Zurich, and Freie Universitaet Berlin 2002-2018.
6 //
7 // This software is released under a three-clause BSD license:
8 // * Redistributions of source code must retain the above copyright
9 // notice, this list of conditions and the following disclaimer.
10 // * Redistributions in binary form must reproduce the above copyright
11 // notice, this list of conditions and the following disclaimer in the
12 // documentation and/or other materials provided with the distribution.
13 // * Neither the name of any author or any participating institution
14 // may be used to endorse or promote products derived from this software
15 // without specific prior written permission.
16 // For a full list of authors, refer to the file AUTHORS.
17 // --------------------------------------------------------------------------
18 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
19 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
21 // ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING
22 // INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
23 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
24 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
25 // OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
26 // WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
27 // OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
28 // ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
29 //
30 // --------------------------------------------------------------------------
31 // $Maintainer: Timo Sachsenberg $
32 // $Authors: Andreas Bertsch $
33 // --------------------------------------------------------------------------
34 //
35 #pragma once
36 
38 
39 #include <vector>
40 #include <map>
41 #include <utility>
42 #include <algorithm>
43 
44 #define ALIGNMENT_DEBUG
45 #undef ALIGNMENT_DEBUG
46 
47 namespace OpenMS
48 {
49 
65  class OPENMS_DLLAPI SpectrumAlignment :
66  public DefaultParamHandler
67  {
68 public:
69 
70  // @name Constructors and Destructors
71  // @{
74 
76  SpectrumAlignment(const SpectrumAlignment & source);
77 
79  ~SpectrumAlignment() override;
80 
82  SpectrumAlignment & operator=(const SpectrumAlignment & source);
83  // @}
84 
85  template <typename SpectrumType1, typename SpectrumType2>
86  void getSpectrumAlignment(std::vector<std::pair<Size, Size> > & alignment, const SpectrumType1 & s1, const SpectrumType2 & s2) const
87  {
88  if (!s1.isSorted() || !s2.isSorted())
89  {
90  throw Exception::IllegalArgument(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "Input to SpectrumAlignment is not sorted!");
91  }
92 
93  // clear result
94  alignment.clear();
95  double tolerance = (double)param_.getValue("tolerance");
96 
97  if (!param_.getValue("is_relative_tolerance").toBool() )
98  {
99  std::map<Size, std::map<Size, std::pair<Size, Size> > > traceback;
100  std::map<Size, std::map<Size, double> > matrix;
101 
102  // init the matrix with "gap costs" tolerance
103  matrix[0][0] = 0;
104  for (Size i = 1; i <= s1.size(); ++i)
105  {
106  matrix[i][0] = i * tolerance;
107  traceback[i][0] = std::make_pair(i - 1, 0);
108  }
109  for (Size j = 1; j <= s2.size(); ++j)
110  {
111  matrix[0][j] = j * tolerance;
112  traceback[0][j] = std::make_pair(0, j - 1);
113  }
114 
115  // fill in the matrix
116  Size left_ptr(1);
117  Size last_i(0), last_j(0);
118 
119  //Size off_band_counter(0);
120  for (Size i = 1; i <= s1.size(); ++i)
121  {
122  double pos1(s1[i - 1].getMZ());
123 
124  for (Size j = left_ptr; j <= s2.size(); ++j)
125  {
126  bool off_band(false);
127  // find min of the three possible directions
128  double pos2(s2[j - 1].getMZ());
129  double diff_align = fabs(pos1 - pos2);
130 
131  // running off the right border of the band?
132  if (pos2 > pos1 && diff_align > tolerance)
133  {
134  if (i < s1.size() && j < s2.size() && s1[i].getMZ() < pos2)
135  {
136  off_band = true;
137  }
138  }
139 
140  // can we tighten the left border of the band?
141  if (pos1 > pos2 && diff_align > tolerance && j > left_ptr + 1)
142  {
143  ++left_ptr;
144  }
145 
146  double score_align = diff_align;
147 
148  if (matrix.find(i - 1) != matrix.end() && matrix[i - 1].find(j - 1) != matrix[i - 1].end())
149  {
150  score_align += matrix[i - 1][j - 1];
151  }
152  else
153  {
154  score_align += (i - 1 + j - 1) * tolerance;
155  }
156 
157  double score_up = tolerance;
158  if (matrix.find(i) != matrix.end() && matrix[i].find(j - 1) != matrix[i].end())
159  {
160  score_up += matrix[i][j - 1];
161  }
162  else
163  {
164  score_up += (i + j - 1) * tolerance;
165  }
166 
167  double score_left = tolerance;
168  if (matrix.find(i - 1) != matrix.end() && matrix[i - 1].find(j) != matrix[i - 1].end())
169  {
170  score_left += matrix[i - 1][j];
171  }
172  else
173  {
174  score_left += (i - 1 + j) * tolerance;
175  }
176 
177  #ifdef ALIGNMENT_DEBUG
178  cerr << i << " " << j << " " << left_ptr << " " << pos1 << " " << pos2 << " " << score_align << " " << score_left << " " << score_up << endl;
179  #endif
180 
181  if (score_align <= score_up && score_align <= score_left && diff_align <= tolerance)
182  {
183  matrix[i][j] = score_align;
184  traceback[i][j] = std::make_pair(i - 1, j - 1);
185  last_i = i;
186  last_j = j;
187  }
188  else
189  {
190  if (score_up <= score_left)
191  {
192  matrix[i][j] = score_up;
193  traceback[i][j] = std::make_pair(i, j - 1);
194  }
195  else
196  {
197  matrix[i][j] = score_left;
198  traceback[i][j] = std::make_pair(i - 1, j);
199  }
200  }
201 
202  if (off_band)
203  {
204  break;
205  }
206  }
207  }
208 
209  //last_i = s1.size() + 1;
210  //last_j = s2.size() + 1;
211 
212  //cerr << last_i << " " << last_j << endl;
213 
214  #ifdef ALIGNMENT_DEBUG
215  #if 0
216  cerr << "TheMatrix: " << endl << " \t \t";
217  for (Size j = 0; j != s2.size(); ++j)
218  {
219  cerr << s2[j].getPosition()[0] << " \t";
220  }
221  cerr << endl;
222  for (Size i = 0; i <= s1.size(); ++i)
223  {
224  if (i != 0)
225  {
226  cerr << s1[i - 1].getPosition()[0] << " \t";
227  }
228  else
229  {
230  cerr << " \t";
231  }
232  for (Size j = 0; j <= s2.size(); ++j)
233  {
234  if (matrix.has(i) && matrix[i].has(j))
235  {
236  if (traceback[i][j].first == i - 1 && traceback[i][j].second == j - 1)
237  {
238  cerr << "\\";
239  }
240  else
241  {
242  if (traceback[i][j].first == i - 1 && traceback[i][j].second == j)
243  {
244  cerr << "|";
245  }
246  else
247  {
248  cerr << "-";
249  }
250  }
251 
252  cerr << matrix[i][j] << " \t";
253  }
254  else
255  {
256  cerr << "-1 \t";
257  }
258  }
259  cerr << endl;
260  }
261  #endif
262  #endif
263 
264  // do traceback
265  Size i = last_i;
266  Size j = last_j;
267 
268  while (i >= 1 && j >= 1)
269  {
270  if (traceback[i][j].first == i - 1 && traceback[i][j].second == j - 1)
271  {
272  alignment.push_back(std::make_pair(i - 1, j - 1));
273  }
274  Size new_i = traceback[i][j].first;
275  Size new_j = traceback[i][j].second;
276 
277  i = new_i;
278  j = new_j;
279  }
280 
281  std::reverse(alignment.begin(), alignment.end());
282 
283  #ifdef ALIGNMENT_DEBUG
284  #if 0
285  // print alignment
286  cerr << "Alignment (size=" << alignment.size() << "): " << endl;
287 
288  Size i_s1(0), i_s2(0);
289  for (vector<pair<Size, Size> >::const_reverse_iterator it = alignment.rbegin(); it != alignment.rend(); ++it, ++i_s1, ++i_s2)
290  {
291  while (i_s1 < it->first - 1)
292  {
293  cerr << i_s1 << " " << s1[i_s1].getPosition()[0] << " " << s1[i_s1].getIntensity() << endl;
294  i_s1++;
295  }
296  while (i_s2 < it->second - 1)
297  {
298  cerr << " \t " << i_s2 << " " << s2[i_s2].getPosition()[0] << " " << s2[i_s2].getIntensity() << endl;
299  i_s2++;
300  }
301  cerr << "(" << s1[it->first - 1].getPosition()[0] << " <-> " << s2[it->second - 1].getPosition()[0] << ") ("
302  << it->first << "|" << it->second << ") (" << s1[it->first - 1].getIntensity() << "|" << s2[it->second - 1].getIntensity() << ")" << endl;
303  }
304  #endif
305  #endif
306  }
307  else // relative alignment (ppm tolerance)
308  {
309  for (Size i = 0; i != s1.size(); ++i)
310  {
311  const double& theo_mz = s1[i].getMZ();
312  double max_dist_dalton = theo_mz * tolerance * 1e-6;
313 
314  // iterate over peaks in experimental spectrum in given fragment tolerance around theoretical peak
315  Size j = s2.findNearest(theo_mz);
316  double exp_mz = s2[j].getMZ();
317 
318  // found peak match
319  if (std::abs(theo_mz - exp_mz) < max_dist_dalton)
320  {
321  alignment.push_back(std::make_pair(i, j));
322  }
323  }
324  }
325  }
326  };
327 }
Main OpenMS namespace.
Definition: FeatureDeconvolution.h:46
A method or algorithm argument contains illegal values.
Definition: Exception.h:648
void getSpectrumAlignment(std::vector< std::pair< Size, Size > > &alignment, const SpectrumType1 &s1, const SpectrumType2 &s2) const
Definition: SpectrumAlignment.h:86
Aligns the peaks of two sorted spectra Method 1: Using a banded (width via &#39;tolerance&#39; parameter) ali...
Definition: SpectrumAlignment.h:65
size_t Size
Size type e.g. used as variable which can hold result of size()
Definition: Types.h:127
A base class for all classes handling default parameters.
Definition: DefaultParamHandler.h:91