OpenMS
SpectrumAlignment.h
Go to the documentation of this file.
1 // Copyright (c) 2002-present, The OpenMS Team -- EKU Tuebingen, ETH Zurich, and FU Berlin
2 // SPDX-License-Identifier: BSD-3-Clause
3 //
4 // --------------------------------------------------------------------------
5 // $Maintainer: Timo Sachsenberg $
6 // $Authors: Andreas Bertsch $
7 // --------------------------------------------------------------------------
8 //
9 #pragma once
10 
13 
14 #include <vector>
15 #include <map>
16 #include <utility>
17 #include <algorithm>
18 
19 #define ALIGNMENT_DEBUG
20 #undef ALIGNMENT_DEBUG
21 
22 namespace OpenMS
23 {
24 
41  class OPENMS_DLLAPI SpectrumAlignment :
42  public DefaultParamHandler
43  {
44 public:
45 
46  // @name Constructors and Destructors
47  // @{
50 
53 
55  ~SpectrumAlignment() override;
56 
59  // @}
60 
61  template <typename SpectrumType1, typename SpectrumType2>
62  void getSpectrumAlignment(std::vector<std::pair<Size, Size> >& alignment, const SpectrumType1& s1, const SpectrumType2& s2) const
63  {
64  if (!s1.isSorted() || !s2.isSorted())
65  {
66  throw Exception::IllegalArgument(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "Input to SpectrumAlignment is not sorted!");
67  }
68 
69  // clear result
70  alignment.clear();
71  double tolerance = (double)param_.getValue("tolerance");
72 
73  if (!param_.getValue("is_relative_tolerance").toBool() )
74  {
75  std::map<Size, std::map<Size, std::pair<Size, Size> > > traceback;
76  std::map<Size, std::map<Size, double> > matrix;
77 
78  // init the matrix with "gap costs" tolerance
79  matrix[0][0] = 0;
80  for (Size i = 1; i <= s1.size(); ++i)
81  {
82  matrix[i][0] = i * tolerance;
83  traceback[i][0] = std::make_pair(i - 1, 0);
84  }
85  for (Size j = 1; j <= s2.size(); ++j)
86  {
87  matrix[0][j] = j * tolerance;
88  traceback[0][j] = std::make_pair(0, j - 1);
89  }
90 
91  // fill in the matrix
92  Size left_ptr(1);
93  Size last_i(0), last_j(0);
94 
95  //Size off_band_counter(0);
96  for (Size i = 1; i <= s1.size(); ++i)
97  {
98  double pos1(s1[i - 1].getMZ());
99 
100  for (Size j = left_ptr; j <= s2.size(); ++j)
101  {
102  bool off_band(false);
103  // find min of the three possible directions
104  double pos2(s2[j - 1].getMZ());
105  double diff_align = fabs(pos1 - pos2);
106 
107  // running off the right border of the band?
108  if (pos2 > pos1 && diff_align > tolerance)
109  {
110  if (i < s1.size() && j < s2.size() && s1[i].getMZ() < pos2)
111  {
112  off_band = true;
113  }
114  }
115 
116  // can we tighten the left border of the band?
117  if (pos1 > pos2 && diff_align > tolerance && j > left_ptr + 1)
118  {
119  ++left_ptr;
120  }
121 
122  double score_align = diff_align;
123 
124  if (matrix.find(i - 1) != matrix.end() && matrix[i - 1].find(j - 1) != matrix[i - 1].end())
125  {
126  score_align += matrix[i - 1][j - 1];
127  }
128  else
129  {
130  score_align += (i - 1 + j - 1) * tolerance;
131  }
132 
133  double score_up = tolerance;
134  if (matrix.find(i) != matrix.end() && matrix[i].find(j - 1) != matrix[i].end())
135  {
136  score_up += matrix[i][j - 1];
137  }
138  else
139  {
140  score_up += (i + j - 1) * tolerance;
141  }
142 
143  double score_left = tolerance;
144  if (matrix.find(i - 1) != matrix.end() && matrix[i - 1].find(j) != matrix[i - 1].end())
145  {
146  score_left += matrix[i - 1][j];
147  }
148  else
149  {
150  score_left += (i - 1 + j) * tolerance;
151  }
152 
153  #ifdef ALIGNMENT_DEBUG
154  cerr << i << " " << j << " " << left_ptr << " " << pos1 << " " << pos2 << " " << score_align << " " << score_left << " " << score_up << endl;
155  #endif
156 
157  if (score_align <= score_up && score_align <= score_left && diff_align <= tolerance)
158  {
159  matrix[i][j] = score_align;
160  traceback[i][j] = std::make_pair(i - 1, j - 1);
161  last_i = i;
162  last_j = j;
163  }
164  else
165  {
166  if (score_up <= score_left)
167  {
168  matrix[i][j] = score_up;
169  traceback[i][j] = std::make_pair(i, j - 1);
170  }
171  else
172  {
173  matrix[i][j] = score_left;
174  traceback[i][j] = std::make_pair(i - 1, j);
175  }
176  }
177 
178  if (off_band)
179  {
180  break;
181  }
182  }
183  }
184 
185  //last_i = s1.size() + 1;
186  //last_j = s2.size() + 1;
187 
188  //cerr << last_i << " " << last_j << endl;
189 
190  #ifdef ALIGNMENT_DEBUG
191  #if 0
192  cerr << "TheMatrix: " << endl << " \t \t";
193  for (Size j = 0; j != s2.size(); ++j)
194  {
195  cerr << s2[j].getPosition()[0] << " \t";
196  }
197  cerr << endl;
198  for (Size i = 0; i <= s1.size(); ++i)
199  {
200  if (i != 0)
201  {
202  cerr << s1[i - 1].getPosition()[0] << " \t";
203  }
204  else
205  {
206  cerr << " \t";
207  }
208  for (Size j = 0; j <= s2.size(); ++j)
209  {
210  if (matrix.has(i) && matrix[i].has(j))
211  {
212  if (traceback[i][j].first == i - 1 && traceback[i][j].second == j - 1)
213  {
214  cerr << "\\";
215  }
216  else
217  {
218  if (traceback[i][j].first == i - 1 && traceback[i][j].second == j)
219  {
220  cerr << "|";
221  }
222  else
223  {
224  cerr << "-";
225  }
226  }
227 
228  cerr << matrix[i][j] << " \t";
229  }
230  else
231  {
232  cerr << "-1 \t";
233  }
234  }
235  cerr << endl;
236  }
237  #endif
238  #endif
239 
240  // do traceback
241  Size i = last_i;
242  Size j = last_j;
243 
244  while (i >= 1 && j >= 1)
245  {
246  if (traceback[i][j].first == i - 1 && traceback[i][j].second == j - 1)
247  {
248  alignment.push_back(std::make_pair(i - 1, j - 1));
249  }
250  Size new_i = traceback[i][j].first;
251  Size new_j = traceback[i][j].second;
252 
253  i = new_i;
254  j = new_j;
255  }
256 
257  std::reverse(alignment.begin(), alignment.end());
258 
259  #ifdef ALIGNMENT_DEBUG
260  #if 0
261  // print alignment
262  cerr << "Alignment (size=" << alignment.size() << "): " << endl;
263 
264  Size i_s1(0), i_s2(0);
265  for (vector<pair<Size, Size> >::const_reverse_iterator it = alignment.rbegin(); it != alignment.rend(); ++it, ++i_s1, ++i_s2)
266  {
267  while (i_s1 < it->first - 1)
268  {
269  cerr << i_s1 << " " << s1[i_s1].getPosition()[0] << " " << s1[i_s1].getIntensity() << endl;
270  i_s1++;
271  }
272  while (i_s2 < it->second - 1)
273  {
274  cerr << " \t " << i_s2 << " " << s2[i_s2].getPosition()[0] << " " << s2[i_s2].getIntensity() << endl;
275  i_s2++;
276  }
277  cerr << "(" << s1[it->first - 1].getPosition()[0] << " <-> " << s2[it->second - 1].getPosition()[0] << ") ("
278  << it->first << "|" << it->second << ") (" << s1[it->first - 1].getIntensity() << "|" << s2[it->second - 1].getIntensity() << ")" << endl;
279  }
280  #endif
281  #endif
282  }
283  else // relative alignment (ppm tolerance)
284  {
285  // find closest match of s1[i] in s2 for all i
286  MatchedIterator<SpectrumType1, PpmTrait> it(s1, s2, tolerance);
287  for (; it != it.end(); ++it) alignment.emplace_back(it.refIdx(), it.tgtIdx());
288  }
289  }
290  };
291 }
A base class for all classes handling default parameters.
Definition: DefaultParamHandler.h:66
A method or algorithm argument contains illegal values.
Definition: Exception.h:616
For each element in the reference container the closest peak in the target will be searched....
Definition: MatchedIterator.h:46
size_t refIdx() const
index into reference container
Definition: MatchedIterator.h:162
size_t tgtIdx() const
index into target container
Definition: MatchedIterator.h:168
static MatchedIterator end()
the end iterator
Definition: MatchedIterator.h:196
Aligns the peaks of two sorted spectra Method 1: Using a banded (width via 'tolerance' parameter) ali...
Definition: SpectrumAlignment.h:43
void getSpectrumAlignment(std::vector< std::pair< Size, Size > > &alignment, const SpectrumType1 &s1, const SpectrumType2 &s2) const
Definition: SpectrumAlignment.h:62
SpectrumAlignment & operator=(const SpectrumAlignment &source)
assignment operator
SpectrumAlignment(const SpectrumAlignment &source)
copy constructor
SpectrumAlignment()
default constructor
~SpectrumAlignment() override
destructor
size_t Size
Size type e.g. used as variable which can hold result of size()
Definition: Types.h:97
static String & reverse(String &this_s)
Definition: StringUtilsSimple.h:330
Main OpenMS namespace.
Definition: openswathalgo/include/OpenMS/OPENSWATHALGO/DATAACCESS/ISpectrumAccess.h:19