libpappsomspp
Library for mass spectrometry
timsframebase.cpp
Go to the documentation of this file.
1 /**
2  * \file pappsomspp/vendors/tims/timsframebase.cpp
3  * \date 16/12/2019
4  * \author Olivier Langella
5  * \brief handle a single Bruker's TimsTof frame without binary data
6  */
7 
8 /*******************************************************************************
9  * Copyright (c) 2019 Olivier Langella <Olivier.Langella@u-psud.fr>.
10  *
11  * This file is part of the PAPPSOms++ library.
12  *
13  * PAPPSOms++ is free software: you can redistribute it and/or modify
14  * it under the terms of the GNU General Public License as published by
15  * the Free Software Foundation, either version 3 of the License, or
16  * (at your option) any later version.
17  *
18  * PAPPSOms++ is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU General Public License for more details.
22  *
23  * You should have received a copy of the GNU General Public License
24  * along with PAPPSOms++. If not, see <http://www.gnu.org/licenses/>.
25  *
26  ******************************************************************************/
27 #include "timsframebase.h"
28 #include "../../../pappsomspp/pappsoexception.h"
29 #include "../../../pappsomspp/exception/exceptionoutofrange.h"
30 #include <QDebug>
31 #include <cmath>
32 #include <solvers.h>
33 
34 namespace pappso
35 {
36 
37 TimsFrameBase::TimsFrameBase(std::size_t timsId, quint32 scanNum)
38 {
39  qDebug() << timsId;
40  m_timsId = timsId;
41 
42  m_scanNumber = scanNum;
43 }
44 
45 TimsFrameBase::TimsFrameBase([[maybe_unused]] const TimsFrameBase &other)
46 {
47 }
48 
50 {
51 }
52 
53 void
54 TimsFrameBase::setAccumulationTime(double accumulation_time_ms)
55 {
56  m_accumulationTime = accumulation_time_ms;
57 }
58 
59 
60 void
61 TimsFrameBase::setMzCalibration(double temperature_correction,
62  double digitizerTimebase,
63  double digitizerDelay,
64  double C0,
65  double C1,
66  double C2,
67  double C3)
68 {
69 
70  // temperature compensation
71  C1 = C1 * temperature_correction;
72  C2 = C2 / temperature_correction;
73 
74  m_digitizerDelay = digitizerDelay;
75  m_digitizerTimebase = digitizerTimebase;
76 
77  m_mzCalibrationArr.push_back(C0);
78  m_mzCalibrationArr.push_back(std::sqrt(std::pow(10, 12) / C1));
79  m_mzCalibrationArr.push_back(C2);
80  m_mzCalibrationArr.push_back(C3);
81 }
82 
83 bool
84 TimsFrameBase::checkScanNum(std::size_t scanNum) const
85 {
86  if(scanNum >= m_scanNumber)
87  {
89  QObject::tr("Invalid scan number : scanNum%1 > m_scanNumber")
90  .arg(scanNum));
91  }
92 
93  return true;
94 }
95 
96 std::size_t
97 TimsFrameBase::getNbrPeaks(std::size_t scanNum) const
98 {
99  throw PappsoException(
100  QObject::tr(
101  "ERROR unable to get number of peaks in TimsFrameBase for scan number %1")
102  .arg(scanNum));
103 }
104 
106 TimsFrameBase::getMassSpectrumSPtr(std::size_t scanNum) const
107 {
108  throw PappsoException(
109  QObject::tr(
110  "ERROR unable to getMassSpectrumSPtr in TimsFrameBase for scan number %1")
111  .arg(scanNum));
112 }
113 Trace
114 TimsFrameBase::cumulateScanToTrace(std::size_t scanNumBegin,
115  std::size_t scanNumEnd) const
116 {
117  throw PappsoException(
118  QObject::tr("ERROR unable to cumulateScanToTrace in TimsFrameBase for scan "
119  "number begin %1 end %2")
120  .arg(scanNumBegin)
121  .arg(scanNumEnd));
122 }
123 void
124 TimsFrameBase::cumulateScansInRawMap(std::map<quint32, quint32> &rawSpectrum
125  [[maybe_unused]],
126  std::size_t scanNumBegin,
127  std::size_t scanNumEnd) const
128 {
129  throw PappsoException(
130  QObject::tr(
131  "ERROR unable to cumulateScansInRawMap in TimsFrameBase for scan "
132  "number begin %1 end %2")
133  .arg(scanNumBegin)
134  .arg(scanNumEnd));
135 }
136 double
138 {
139  // mz calibration
140  return (index * m_digitizerTimebase) + m_digitizerDelay;
141 }
142 double
143 TimsFrameBase::getTofFromIndex(quint32 index) const
144 {
145  // mz calibration
146  return ((double)index * m_digitizerTimebase) + m_digitizerDelay;
147 }
148 double
150 {
151  // http://www.alglib.net/equations/polynomial.php
152  // http://www.alglib.net/translator/man/manual.cpp.html#sub_polynomialsolve
153  // https://math.stackexchange.com/questions/1291208/number-of-roots-of-a-polynomial-of-non-integer-degree
154  // https://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=2&ved=2ahUKEwiWhLOFxqrkAhVLxYUKHVqqDFcQFjABegQIAxAB&url=https%3A%2F%2Fkluge.in-chemnitz.de%2Fopensource%2Fspline%2Fexample_alglib.cpp&usg=AOvVaw0guGejJGPmkOVg48_GJYR8
155  // https://stackoverflow.com/questions/26091323/how-to-plot-a-function-curve-in-r
156  /*
157  * beware to put the function on a single line in R:
158 > eq <- function(m){ 1 + (sqrt((10^12)/670) * sqrt(m)) + (207.775676931964 * m)
159 + (59.2526676368822 * (m^1.5)) }
160 > eq <- function(m){ 313.577620892277 + (sqrt((10^12)/157424.07710945) *
161 sqrt(m)) + (0.000338743021989553 * m)
162 + (0 * (m^1.5)) }
163 > plot(eq(1:1000), type='l')
164 
165 
166 
167 > eq2 <- function(m2){ 1 + sqrt((10^12)/670) * m2 + 207.775676931964 * (m2^2)
168 + 59.2526676368822 * (m2^3) }
169 > plot(eq2(1:sqrt(1000)), type='l')
170 */
171  // How to Factor a Trinomial with Fractions as Coefficients
172 
173  // formula
174  // a = c0 = 1
175  // b = sqrt((10^12)/c1), c1 = 670 * m^0.5 (1/2)
176  // c = c2, c2 = 207.775676931964 * m
177  // d = c3, c3 = 59.2526676368822 * m^1.5 (3/2)
178  // double mz = 0;
179  std::vector<double> X;
180  X.push_back((m_mzCalibrationArr[0] - (double)tof));
181  X.push_back(m_mzCalibrationArr[1]);
182  if(m_mzCalibrationArr[2] != 0)
183  {
184  X.push_back(m_mzCalibrationArr[2]);
185  }
186  if(m_mzCalibrationArr[3] != 0)
187  {
188  X.push_back(m_mzCalibrationArr[3]);
189  }
190  else
191  {
192  qDebug() << m_mzCalibrationArr[3];
193  }
194 
195  alglib::real_1d_array polynom_array;
196  try
197  {
198  polynom_array.setcontent(X.size(), &(X[0]));
199  }
200  catch(alglib::ap_error &error)
201  {
202  // PolynomialSolve: A[N]=0
204  QObject::tr("ERROR in alglib::polynom_array.setcontent :\n%1")
205  .arg(error.msg.c_str()));
206  }
207 
208 
209  /*
210  alglib::polynomialsolve(
211 real_1d_array a,
212 ae_int_t n,
213 complex_1d_array& x,
214 polynomialsolverreport& rep,
215 const xparams _params = alglib::xdefault);
216 */
217  alglib::complex_1d_array m;
218  alglib::polynomialsolverreport rep;
219  // qDebug();
220  try
221  {
222  alglib::polynomialsolve(polynom_array, X.size() - 1, m, rep);
223  }
224  catch(alglib::ap_error &error)
225  {
226  qDebug() << " X.size() - 1 = " << X.size() - 1;
227  qDebug() << m_mzCalibrationArr[0];
228  qDebug() << m_mzCalibrationArr[1];
229  qDebug() << m_mzCalibrationArr[2];
230  qDebug() << m_mzCalibrationArr[3];
231 
232  // PolynomialSolve: A[N]=0
234  QObject::tr("ERROR in alglib::polynomialsolve :\n%1")
235  .arg(error.msg.c_str()));
236  }
237 
238 
239  // qDebug();
240 
241  if(m.length() == 0)
242  {
244  QObject::tr("ERROR in TimsFrame::getMzFromTof m.size() == 0"));
245  }
246  // qDebug();
247  if(m[0].y != 0)
248  {
250  QObject::tr("ERROR in TimsFrame::getMzFromTof m[0].y!= 0"));
251  }
252 
253  return pow(m[0].x, 2);
254 }
255 
256 quint32
258 {
259  // formula
260  // a = c0 = 1
261  // b = sqrt((10^12)/c1), c1 = 670 * m^0.5 (1/2)
262  // c = c2, c2 = 207.775676931964 * m
263  // d = c3, c3 = 59.2526676368822 * m^1.5 (3/2)
264  qDebug() << "mz=" << mz;
265 
266  double tof = m_mzCalibrationArr[0];
267  qDebug() << "tof ( m_mzCalibrationArr[0])=" << tof;
268  // TODO cache value of std::sqrt((std::pow(10, 12) / m_mzCalibrationArr[1]))
269  tof += m_mzCalibrationArr[1] * std::sqrt(mz);
270  qDebug() << "tof=" << tof;
271  tof += m_mzCalibrationArr[2] * mz;
272  qDebug() << "tof=" << tof;
273  tof += m_mzCalibrationArr[3] * std::pow(mz, 1.5);
274  qDebug() << "tof=" << tof;
275  tof -= m_digitizerDelay;
276  qDebug() << "tof=" << tof;
277  tof = tof / m_digitizerTimebase;
278  qDebug() << "index=" << tof;
279  return (quint32)std::round(tof);
280 }
281 
282 void
284 {
285  m_time = time;
286 }
287 
288 void
290 {
291 
292  qDebug() << " m_msMsType=" << type;
293  m_msMsType = type;
294 }
295 
296 unsigned int
298 {
299  if(m_msMsType == 0)
300  return 1;
301  return 2;
302 }
303 
304 double
306 {
307  return m_time;
308 }
309 
310 std::size_t
312 {
313  return m_timsId;
314 }
315 void
317  double C0,
318  double C1,
319  double C2,
320  double C3,
321  double C4,
322  [[maybe_unused]] double C5,
323  double C6,
324  double C7,
325  double C8,
326  double C9)
327 {
328  if(tims_model_type != 2)
329  {
330  throw pappso::PappsoException(QObject::tr(
331  "ERROR in TimsFrame::setTimsCalibration tims_model_type != 2"));
332  }
333  m_timsDvStart = C2; // C2 from TimsCalibration
334  m_timsTtrans = C4; // C4 from TimsCalibration
335  m_timsNdelay = C0; // C0 from TimsCalibration
336  m_timsVmin = C8; // C8 from TimsCalibration
337  m_timsVmax = C9; // C9 from TimsCalibration
338  m_timsC6 = C6;
339  m_timsC7 = C7;
340 
341 
342  m_timsSlope =
343  (C3 - m_timsDvStart) / C1; // //C3 from TimsCalibration // C2 from
344  // TimsCalibration // C1 from TimsCalibration
345 }
346 double
347 TimsFrameBase::getVoltageTransformation(std::size_t scanNum) const
348 {
349  double v = m_timsDvStart +
350  m_timsSlope * ((double)scanNum - m_timsTtrans - m_timsNdelay);
351 
352  if(v < m_timsVmin)
353  {
355  QObject::tr("ERROR in TimsFrame::getVoltageTransformation invalid tims "
356  "calibration, v < m_timsVmin"));
357  }
358 
359 
360  if(v > m_timsVmax)
361  {
363  QObject::tr("ERROR in TimsFrame::getVoltageTransformation invalid tims "
364  "calibration, v > m_timsVmax"));
365  }
366  return v;
367 }
368 double
369 TimsFrameBase::getDriftTime(std::size_t scanNum) const
370 {
371  return (m_accumulationTime / (double)m_scanNumber) * ((double)scanNum);
372 }
373 
374 double
376 {
377  return 1 / (m_timsC6 + (m_timsC7 / getVoltageTransformation(scanNum)));
378 }
379 
380 
381 std::size_t
382 TimsFrameBase::getScanNumFromOneOverK0(double one_over_k0) const
383 {
384  double temp = 1 / one_over_k0;
385  temp = temp - m_timsC6;
386  temp = temp / m_timsC7;
387  temp = 1 / temp;
388  temp = temp - m_timsDvStart;
389  temp = temp / m_timsSlope + m_timsTtrans + m_timsNdelay;
390  return (std::size_t)std::round(temp);
391 }
392 
393 bool
395 {
396  if((m_timsDvStart == other.m_timsDvStart) &&
397  (m_timsTtrans == other.m_timsTtrans) &&
398  (m_timsNdelay == other.m_timsNdelay) && (m_timsVmin == other.m_timsVmin) &&
399  (m_timsVmax == other.m_timsVmax) && (m_timsC6 == other.m_timsC6) &&
400  (m_timsC7 == other.m_timsC7) && (m_timsSlope == other.m_timsSlope))
401  {
402  return true;
403  }
404  return false;
405 }
406 
407 
410  std::map<quint32, quint32> &accumulated_scans) const
411 {
412  qDebug();
413  // qDebug();
414  // add flanking peaks
415  pappso::Trace local_trace;
416  transform(begin(accumulated_scans),
417  end(accumulated_scans),
418  back_inserter(local_trace),
419  [](std::map<quint32, quint32>::value_type const &pair) {
420  return pappso::DataPoint(pair.first, pair.second);
421  });
422 
423 
424  for(DataPoint &element : local_trace)
425  {
426  // intensity normalization
427  element.y *= 100.0 / m_accumulationTime;
428 
429  // mz calibration
430  double tof = (element.x * m_digitizerTimebase) + m_digitizerDelay;
431  element.x = getMzFromTof(tof);
432  }
433  local_trace.sortX();
434 
435  qDebug();
436  // qDebug();
437  return local_trace;
438 }
439 
442  std::map<quint32, quint32> &accumulated_scans) const
443 {
444  qDebug();
445  // qDebug();
446  // add flanking peaks
447  std::vector<quint32> keys;
448  transform(begin(accumulated_scans),
449  end(accumulated_scans),
450  back_inserter(keys),
451  [](std::map<quint32, quint32>::value_type const &pair) {
452  return pair.first;
453  });
454  std::sort(keys.begin(), keys.end());
455  pappso::DataPoint data_point_cumul;
456  data_point_cumul.x = 0;
457  data_point_cumul.y = 0;
458 
459  pappso::Trace local_trace;
460 
461  quint32 last_key = 0;
462 
463  for(quint32 key : keys)
464  {
465  if(key == last_key + 1)
466  {
467  // cumulate
468  if(accumulated_scans[key] > accumulated_scans[last_key])
469  {
470  if(data_point_cumul.x == last_key)
471  {
472  // growing peak
473  data_point_cumul.x = key;
474  data_point_cumul.y += accumulated_scans[key];
475  }
476  else
477  {
478  // new peak
479  // flush
480  if(data_point_cumul.y > 0)
481  {
482  // intensity normalization
483  data_point_cumul.y *= 100.0 / m_accumulationTime;
484 
485 
486  // mz calibration
487  double tof = (data_point_cumul.x * m_digitizerTimebase) +
489  data_point_cumul.x = getMzFromTof(tof);
490  local_trace.push_back(data_point_cumul);
491  }
492 
493  // new point
494  data_point_cumul.x = key;
495  data_point_cumul.y = accumulated_scans[key];
496  }
497  }
498  else
499  {
500  data_point_cumul.y += accumulated_scans[key];
501  }
502  }
503  else
504  {
505  // flush
506  if(data_point_cumul.y > 0)
507  {
508  // intensity normalization
509  data_point_cumul.y *= 100.0 / m_accumulationTime;
510 
511 
512  // mz calibration
513  double tof =
514  (data_point_cumul.x * m_digitizerTimebase) + m_digitizerDelay;
515  data_point_cumul.x = getMzFromTof(tof);
516  local_trace.push_back(data_point_cumul);
517  }
518 
519  // new point
520  data_point_cumul.x = key;
521  data_point_cumul.y = accumulated_scans[key];
522  }
523 
524  last_key = key;
525  }
526  // flush
527  if(data_point_cumul.y > 0)
528  {
529  // intensity normalization
530  data_point_cumul.y *= 100.0 / m_accumulationTime;
531 
532 
533  // mz calibration
534  double tof =
535  (data_point_cumul.x * m_digitizerTimebase) + m_digitizerDelay;
536  data_point_cumul.x = getMzFromTof(tof);
537  local_trace.push_back(data_point_cumul);
538  }
539 
540  local_trace.sortX();
541  qDebug();
542  // qDebug();
543  return local_trace;
544 }
545 
546 } // namespace pappso
pappso::TimsFrameBase::m_digitizerTimebase
double m_digitizerTimebase
Definition: timsframebase.h:185
pappso::TimsFrameBase::checkScanNum
bool checkScanNum(std::size_t scanNum) const
Definition: timsframebase.cpp:84
pappso::TimsFrameBase::m_timsId
std::size_t m_timsId
Tims frame database id (the SQL identifier of this frame)
Definition: timsframebase.h:179
pappso::TimsFrameBase::setTimsCalibration
void setTimsCalibration(int tims_model_type, double C0, double C1, double C2, double C3, double C4, double C5, double C6, double C7, double C8, double C9)
Definition: timsframebase.cpp:316
pappso::TimsFrameBase::getMsLevel
unsigned int getMsLevel() const
Definition: timsframebase.cpp:297
pappso::TimsFrameBase::cumulateScanToTrace
virtual Trace cumulateScanToTrace(std::size_t scanNumBegin, std::size_t scanNumEnd) const
Definition: timsframebase.cpp:114
pappso::TimsFrameBase
Definition: timsframebase.h:47
pappso::DataPoint::y
pappso_double y
Definition: datapoint.h:23
pappso::TimsFrameBase::m_timsSlope
double m_timsSlope
Definition: timsframebase.h:199
pappso::TimsFrameBase::m_timsC7
double m_timsC7
Definition: timsframebase.h:207
pappso
tries to keep as much as possible monoisotopes, removing any possible C13 peaks and changes multichar...
Definition: aa.cpp:39
pappso::TimsFrameBase::setAccumulationTime
void setAccumulationTime(double accumulation_time_ms)
Definition: timsframebase.cpp:54
pappso::TimsFrameBase::setMzCalibration
void setMzCalibration(double temperature_correction, double digitizerTimebase, double digitizerDelay, double C0, double C1, double C2, double C3)
Definition: timsframebase.cpp:61
pappso::DataPoint
Definition: datapoint.h:21
pappso::PeptideIonCter::y
@ y
pappso::TimsFrameBase::setTime
void setTime(double time)
Definition: timsframebase.cpp:283
pappso::TimsFrameBase::getNbrPeaks
virtual std::size_t getNbrPeaks(std::size_t scanNum) const
Definition: timsframebase.cpp:97
pappso::TimsFrameBase::getRawIndexFromMz
quint32 getRawIndexFromMz(double mz) const
get raw index of a given m/z
Definition: timsframebase.cpp:257
pappso::TimsFrameBase::getDriftTime
double getDriftTime(std::size_t scanNum) const
get drift time of a scan number in milliseconds
Definition: timsframebase.cpp:369
pappso::TimsFrameBase::m_accumulationTime
double m_accumulationTime
accumulation time in milliseconds
Definition: timsframebase.h:183
pappso::TimsFrameBase::m_time
double m_time
retention time
Definition: timsframebase.h:196
pappso::TimsFrameBase::getTraceFromCumulatedScansBuiltinCentroid
pappso::Trace getTraceFromCumulatedScansBuiltinCentroid(std::map< quint32, quint32 > &accumulated_scans) const
transform accumulation of raw scans into a real mass spectrum with a simple centroid on raw integers
Definition: timsframebase.cpp:441
pappso::ExceptionOutOfRange
Definition: exceptionoutofrange.h:32
pappso::Trace
A simple container of DataPoint instances.
Definition: trace.h:132
pappso::TimsFrameBase::getVoltageTransformation
double getVoltageTransformation(std::size_t scanNum) const
Definition: timsframebase.cpp:347
pappso::TimsFrameBase::getId
std::size_t getId() const
Definition: timsframebase.cpp:311
pappso::TimsFrameBase::m_timsDvStart
double m_timsDvStart
Definition: timsframebase.h:198
pappso::TimsFrameBase::getMassSpectrumSPtr
virtual MassSpectrumSPtr getMassSpectrumSPtr(std::size_t scanNum) const
Definition: timsframebase.cpp:106
pappso::TimsFrameBase::getScanNumFromOneOverK0
std::size_t getScanNumFromOneOverK0(double one_over_k0) const
get the scan number from a given 1/Ko mobility value
Definition: timsframebase.cpp:382
timsframebase.h
handle a single Bruker's TimsTof frame without binary data
pappso::DataPoint::x
pappso_double x
Definition: datapoint.h:22
pappso::Trace::sortX
void sortX()
Definition: trace.cpp:790
pappso::TimsFrameBase::~TimsFrameBase
~TimsFrameBase()
Definition: timsframebase.cpp:49
pappso::TimsFrameBase::TimsFrameBase
TimsFrameBase(std::size_t timsId, quint32 scanNum)
constructor for binary independant tims frame
Definition: timsframebase.cpp:37
pappso::TimsFrameBase::setMsMsType
void setMsMsType(quint8 type)
Definition: timsframebase.cpp:289
pappso::TimsFrameBase::m_timsNdelay
double m_timsNdelay
Definition: timsframebase.h:203
pappso::TimsFrameBase::m_timsVmin
double m_timsVmin
Definition: timsframebase.h:204
pappso::TimsFrameBase::cumulateScansInRawMap
virtual void cumulateScansInRawMap(std::map< quint32, quint32 > &rawSpectrum, std::size_t scanNumBegin, std::size_t scanNumEnd) const
cumulate scan list into a trace into a raw spectrum map
Definition: timsframebase.cpp:124
pappso::TimsFrameBase::getOneOverK0Transformation
double getOneOverK0Transformation(std::size_t scanNum) const
get 1/K0 value of a given scan (mobility value)
Definition: timsframebase.cpp:375
pappso::DataKind::unset
@ unset
not set
pappso::TimsFrameBase::getMzFromTof
double getMzFromTof(double tof) const
get m/z from time of flight
Definition: timsframebase.cpp:149
pappso::TimsFrameBase::m_scanNumber
quint32 m_scanNumber
total number of scans contained in this frame
Definition: timsframebase.h:173
pappso::TimsFrameBase::m_timsC6
double m_timsC6
Definition: timsframebase.h:206
pappso::TimsFrameBase::getTofFromIndex
double getTofFromIndex(quint32 index) const
get time of flight from raw index
Definition: timsframebase.cpp:143
pappso::TimsFrameBase::m_digitizerDelay
double m_digitizerDelay
Definition: timsframebase.h:186
pappso::TimsFrameBase::m_mzCalibrationArr
std::vector< double > m_mzCalibrationArr
MZ calibration parameters.
Definition: timsframebase.h:190
pappso::TimsFrameBase::m_msMsType
quint8 m_msMsType
Definition: timsframebase.h:192
pappso::TimsFrameBase::m_timsVmax
double m_timsVmax
Definition: timsframebase.h:205
pappso::TimsFrameBase::m_timsTtrans
double m_timsTtrans
Definition: timsframebase.h:202
pappso::TimsFrameBase::getTime
double getTime() const
Definition: timsframebase.cpp:305
pappso::MassSpectrumSPtr
std::shared_ptr< MassSpectrum > MassSpectrumSPtr
Definition: massspectrum.h:54
pappso::TimsFrameBase::getTraceFromCumulatedScans
pappso::Trace getTraceFromCumulatedScans(std::map< quint32, quint32 > &accumulated_scans) const
transform accumulation of raw scans into a real mass spectrum
Definition: timsframebase.cpp:409
pappso::TimsFrameBase::hasSameCalibrationData
virtual bool hasSameCalibrationData(const TimsFrameBase &other) const
tells if 2 tims frame has the same calibration data Usefull to know if raw data can be handled betwee...
Definition: timsframebase.cpp:394
pappso::PappsoException
Definition: pappsoexception.h:42