OpenMS  3.0.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-2022.
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 
39 
40 #include <vector>
41 #include <map>
42 #include <utility>
43 #include <algorithm>
44 
45 #define ALIGNMENT_DEBUG
46 #undef ALIGNMENT_DEBUG
47 
48 namespace OpenMS
49 {
50 
67  class OPENMS_DLLAPI SpectrumAlignment :
68  public DefaultParamHandler
69  {
70 public:
71 
72  // @name Constructors and Destructors
73  // @{
76 
78  SpectrumAlignment(const SpectrumAlignment & source);
79 
81  ~SpectrumAlignment() override;
82 
84  SpectrumAlignment & operator=(const SpectrumAlignment & source);
85  // @}
86 
87  template <typename SpectrumType1, typename SpectrumType2>
88  void getSpectrumAlignment(std::vector<std::pair<Size, Size> >& alignment, const SpectrumType1& s1, const SpectrumType2& s2) const
89  {
90  if (!s1.isSorted() || !s2.isSorted())
91  {
92  throw Exception::IllegalArgument(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "Input to SpectrumAlignment is not sorted!");
93  }
94 
95  // clear result
96  alignment.clear();
97  double tolerance = (double)param_.getValue("tolerance");
98 
99  if (!param_.getValue("is_relative_tolerance").toBool() )
100  {
101  std::map<Size, std::map<Size, std::pair<Size, Size> > > traceback;
102  std::map<Size, std::map<Size, double> > matrix;
103 
104  // init the matrix with "gap costs" tolerance
105  matrix[0][0] = 0;
106  for (Size i = 1; i <= s1.size(); ++i)
107  {
108  matrix[i][0] = i * tolerance;
109  traceback[i][0] = std::make_pair(i - 1, 0);
110  }
111  for (Size j = 1; j <= s2.size(); ++j)
112  {
113  matrix[0][j] = j * tolerance;
114  traceback[0][j] = std::make_pair(0, j - 1);
115  }
116 
117  // fill in the matrix
118  Size left_ptr(1);
119  Size last_i(0), last_j(0);
120 
121  //Size off_band_counter(0);
122  for (Size i = 1; i <= s1.size(); ++i)
123  {
124  double pos1(s1[i - 1].getMZ());
125 
126  for (Size j = left_ptr; j <= s2.size(); ++j)
127  {
128  bool off_band(false);
129  // find min of the three possible directions
130  double pos2(s2[j - 1].getMZ());
131  double diff_align = fabs(pos1 - pos2);
132 
133  // running off the right border of the band?
134  if (pos2 > pos1 && diff_align > tolerance)
135  {
136  if (i < s1.size() && j < s2.size() && s1[i].getMZ() < pos2)
137  {
138  off_band = true;
139  }
140  }
141 
142  // can we tighten the left border of the band?
143  if (pos1 > pos2 && diff_align > tolerance && j > left_ptr + 1)
144  {
145  ++left_ptr;
146  }
147 
148  double score_align = diff_align;
149 
150  if (matrix.find(i - 1) != matrix.end() && matrix[i - 1].find(j - 1) != matrix[i - 1].end())
151  {
152  score_align += matrix[i - 1][j - 1];
153  }
154  else
155  {
156  score_align += (i - 1 + j - 1) * tolerance;
157  }
158 
159  double score_up = tolerance;
160  if (matrix.find(i) != matrix.end() && matrix[i].find(j - 1) != matrix[i].end())
161  {
162  score_up += matrix[i][j - 1];
163  }
164  else
165  {
166  score_up += (i + j - 1) * tolerance;
167  }
168 
169  double score_left = tolerance;
170  if (matrix.find(i - 1) != matrix.end() && matrix[i - 1].find(j) != matrix[i - 1].end())
171  {
172  score_left += matrix[i - 1][j];
173  }
174  else
175  {
176  score_left += (i - 1 + j) * tolerance;
177  }
178 
179  #ifdef ALIGNMENT_DEBUG
180  cerr << i << " " << j << " " << left_ptr << " " << pos1 << " " << pos2 << " " << score_align << " " << score_left << " " << score_up << endl;
181  #endif
182 
183  if (score_align <= score_up && score_align <= score_left && diff_align <= tolerance)
184  {
185  matrix[i][j] = score_align;
186  traceback[i][j] = std::make_pair(i - 1, j - 1);
187  last_i = i;
188  last_j = j;
189  }
190  else
191  {
192  if (score_up <= score_left)
193  {
194  matrix[i][j] = score_up;
195  traceback[i][j] = std::make_pair(i, j - 1);
196  }
197  else
198  {
199  matrix[i][j] = score_left;
200  traceback[i][j] = std::make_pair(i - 1, j);
201  }
202  }
203 
204  if (off_band)
205  {
206  break;
207  }
208  }
209  }
210 
211  //last_i = s1.size() + 1;
212  //last_j = s2.size() + 1;
213 
214  //cerr << last_i << " " << last_j << endl;
215 
216  #ifdef ALIGNMENT_DEBUG
217  #if 0
218  cerr << "TheMatrix: " << endl << " \t \t";
219  for (Size j = 0; j != s2.size(); ++j)
220  {
221  cerr << s2[j].getPosition()[0] << " \t";
222  }
223  cerr << endl;
224  for (Size i = 0; i <= s1.size(); ++i)
225  {
226  if (i != 0)
227  {
228  cerr << s1[i - 1].getPosition()[0] << " \t";
229  }
230  else
231  {
232  cerr << " \t";
233  }
234  for (Size j = 0; j <= s2.size(); ++j)
235  {
236  if (matrix.has(i) && matrix[i].has(j))
237  {
238  if (traceback[i][j].first == i - 1 && traceback[i][j].second == j - 1)
239  {
240  cerr << "\\";
241  }
242  else
243  {
244  if (traceback[i][j].first == i - 1 && traceback[i][j].second == j)
245  {
246  cerr << "|";
247  }
248  else
249  {
250  cerr << "-";
251  }
252  }
253 
254  cerr << matrix[i][j] << " \t";
255  }
256  else
257  {
258  cerr << "-1 \t";
259  }
260  }
261  cerr << endl;
262  }
263  #endif
264  #endif
265 
266  // do traceback
267  Size i = last_i;
268  Size j = last_j;
269 
270  while (i >= 1 && j >= 1)
271  {
272  if (traceback[i][j].first == i - 1 && traceback[i][j].second == j - 1)
273  {
274  alignment.push_back(std::make_pair(i - 1, j - 1));
275  }
276  Size new_i = traceback[i][j].first;
277  Size new_j = traceback[i][j].second;
278 
279  i = new_i;
280  j = new_j;
281  }
282 
283  std::reverse(alignment.begin(), alignment.end());
284 
285  #ifdef ALIGNMENT_DEBUG
286  #if 0
287  // print alignment
288  cerr << "Alignment (size=" << alignment.size() << "): " << endl;
289 
290  Size i_s1(0), i_s2(0);
291  for (vector<pair<Size, Size> >::const_reverse_iterator it = alignment.rbegin(); it != alignment.rend(); ++it, ++i_s1, ++i_s2)
292  {
293  while (i_s1 < it->first - 1)
294  {
295  cerr << i_s1 << " " << s1[i_s1].getPosition()[0] << " " << s1[i_s1].getIntensity() << endl;
296  i_s1++;
297  }
298  while (i_s2 < it->second - 1)
299  {
300  cerr << " \t " << i_s2 << " " << s2[i_s2].getPosition()[0] << " " << s2[i_s2].getIntensity() << endl;
301  i_s2++;
302  }
303  cerr << "(" << s1[it->first - 1].getPosition()[0] << " <-> " << s2[it->second - 1].getPosition()[0] << ") ("
304  << it->first << "|" << it->second << ") (" << s1[it->first - 1].getIntensity() << "|" << s2[it->second - 1].getIntensity() << ")" << endl;
305  }
306  #endif
307  #endif
308  }
309  else // relative alignment (ppm tolerance)
310  {
311  // find closest match of s1[i] in s2 for all i
312  MatchedIterator<SpectrumType1, PpmTrait> it(s1, s2, tolerance);
313  for (; it != it.end(); ++it) alignment.emplace_back(it.refIdx(), it.tgtIdx());
314  }
315  }
316  };
317 }
DefaultParamHandler.h
OpenMS::SpectrumAlignment::getSpectrumAlignment
void getSpectrumAlignment(std::vector< std::pair< Size, Size > > &alignment, const SpectrumType1 &s1, const SpectrumType2 &s2) const
Definition: SpectrumAlignment.h:88
double
OpenMS::Exception::IllegalArgument
A method or algorithm argument contains illegal values.
Definition: Exception.h:648
OpenMS::StringUtils::reverse
static String & reverse(String &this_s)
Definition: StringUtilsSimple.h:350
OpenMS::Size
size_t Size
Size type e.g. used as variable which can hold result of size()
Definition: Types.h:127
OpenMS::MatchedIterator
For each element in the reference container the closest peak in the target will be searched....
Definition: MatchedIterator.h:71
MatchedIterator.h
OpenMS::DefaultParamHandler
A base class for all classes handling default parameters.
Definition: DefaultParamHandler.h:92
OpenMS::SpectrumAlignment
Aligns the peaks of two sorted spectra Method 1: Using a banded (width via 'tolerance' parameter) ali...
Definition: SpectrumAlignment.h:67
OpenMS
Main OpenMS namespace.
Definition: FeatureDeconvolution.h:47
OpenMS::MatchedIterator::tgtIdx
size_t tgtIdx() const
index into target container
Definition: MatchedIterator.h:194
OpenMS::MatchedIterator::refIdx
size_t refIdx() const
index into reference container
Definition: MatchedIterator.h:188
OpenMS::MatchedIterator::end
static MatchedIterator end()
the end iterator
Definition: MatchedIterator.h:222