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"
31 #include <QDebug>
32 #include <QObject>
33 #include <cmath>
34 #include <algorithm>
35 
36 namespace pappso
37 {
38 
39 TimsFrameBase::TimsFrameBase(std::size_t timsId, quint32 scanNum)
40 {
41  qDebug() << timsId;
42  m_timsId = timsId;
43 
44  m_scanNumber = scanNum;
45 }
46 
47 TimsFrameBase::TimsFrameBase([[maybe_unused]] const TimsFrameBase &other)
48 {
49 }
50 
52 {
53 }
54 
55 void
56 TimsFrameBase::setAccumulationTime(double accumulation_time_ms)
57 {
58  m_accumulationTime = accumulation_time_ms;
59 }
60 
61 
62 void
64  double T2_frame,
65  double digitizerTimebase,
66  double digitizerDelay,
67  double C0,
68  double C1,
69  double C2,
70  double C3,
71  double C4,
72  double T1_ref,
73  double T2_ref,
74  double dC1,
75  double dC2)
76 {
77 
78  /* MzCalibrationModel1 mzCalibration(temperature_correction,
79  digitizerTimebase,
80  digitizerDelay,
81  C0,
82  C1,
83  C2,
84  C3,
85  C4);
86  */
87  msp_mzCalibration = std::make_shared<MzCalibrationModel1>(T1_frame,
88  T2_frame,
89  digitizerTimebase,
90  digitizerDelay,
91  C0,
92  C1,
93  C2,
94  C3,
95  C4,
96  T1_ref,
97  T2_ref,
98  dC1,
99  dC2);
100 }
101 
102 bool
103 TimsFrameBase::checkScanNum(std::size_t scanNum) const
104 {
105  if(scanNum >= m_scanNumber)
106  {
108  QObject::tr("Invalid scan number : scanNum %1 > m_scanNumber %2")
109  .arg(scanNum)
110  .arg(m_scanNumber));
111  }
112 
113  return true;
114 }
115 
116 std::size_t
117 TimsFrameBase::getNbrPeaks(std::size_t scanNum) const
118 {
119  throw PappsoException(
120  QObject::tr(
121  "ERROR unable to get number of peaks in TimsFrameBase for scan number %1")
122  .arg(scanNum));
123 }
124 
125 std::size_t
127 {
128  return m_scanNumber;
129 }
130 
132 TimsFrameBase::getMassSpectrumSPtr(std::size_t scanNum) const
133 {
134  throw PappsoException(
135  QObject::tr(
136  "ERROR unable to getMassSpectrumSPtr in TimsFrameBase for scan number %1")
137  .arg(scanNum));
138 }
139 Trace
140 TimsFrameBase::cumulateScanToTrace(std::size_t scanNumBegin,
141  std::size_t scanNumEnd) const
142 {
143  throw PappsoException(
144  QObject::tr("ERROR unable to cumulateScanToTrace in TimsFrameBase for scan "
145  "number begin %1 end %2")
146  .arg(scanNumBegin)
147  .arg(scanNumEnd));
148 }
149 
150 void
151 TimsFrameBase::cumulateScansInRawMap(std::map<quint32, quint32> &rawSpectrum
152  [[maybe_unused]],
153  std::size_t scanNumBegin,
154  std::size_t scanNumEnd) const
155 {
156  throw PappsoException(
157  QObject::tr(
158  "ERROR unable to cumulateScansInRawMap in TimsFrameBase for scan "
159  "number begin %1 end %2")
160  .arg(scanNumBegin)
161  .arg(scanNumEnd));
162 }
163 
164 
165 quint64
167 {
168  throw PappsoException(
169  QObject::tr(
170  "ERROR unable to cumulateSingleScanIntensities in TimsFrameBase for scan "
171  "number %1.").arg(scanNum));
172 
173  return 0;
174 }
175 
176 
177 quint64
179  std::size_t scanNumEnd) const
180 {
181  throw PappsoException(
182  QObject::tr(
183  "ERROR unable to cumulateScansInRawMap in TimsFrameBase for scan "
184  "number begin %1 end %2")
185  .arg(scanNumBegin)
186  .arg(scanNumEnd));
187 
188  return 0;
189 }
190 
191 void
193 {
194  m_time = time;
195 }
196 
197 void
199 {
200 
201  qDebug() << " m_msMsType=" << type;
202  m_msMsType = type;
203 }
204 
205 unsigned int
207 {
208  if(m_msMsType == 0)
209  return 1;
210  return 2;
211 }
212 
213 double
215 {
216  return m_time;
217 }
218 
219 std::size_t
221 {
222  return m_timsId;
223 }
224 void
226  double C0,
227  double C1,
228  double C2,
229  double C3,
230  double C4,
231  [[maybe_unused]] double C5,
232  double C6,
233  double C7,
234  double C8,
235  double C9)
236 {
237  if(tims_model_type != 2)
238  {
239  throw pappso::PappsoException(QObject::tr(
240  "ERROR in TimsFrame::setTimsCalibration tims_model_type != 2"));
241  }
242  m_timsDvStart = C2; // C2 from TimsCalibration
243  m_timsTtrans = C4; // C4 from TimsCalibration
244  m_timsNdelay = C0; // C0 from TimsCalibration
245  m_timsVmin = C8; // C8 from TimsCalibration
246  m_timsVmax = C9; // C9 from TimsCalibration
247  m_timsC6 = C6;
248  m_timsC7 = C7;
249 
250 
251  m_timsSlope =
252  (C3 - m_timsDvStart) / C1; // //C3 from TimsCalibration // C2 from
253  // TimsCalibration // C1 from TimsCalibration
254 }
255 double
256 TimsFrameBase::getVoltageTransformation(std::size_t scanNum) const
257 {
258  double v = m_timsDvStart +
259  m_timsSlope * ((double)scanNum - m_timsTtrans - m_timsNdelay);
260 
261  if(v < m_timsVmin)
262  {
264  QObject::tr("ERROR in TimsFrame::getVoltageTransformation invalid tims "
265  "calibration, v < m_timsVmin"));
266  }
267 
268 
269  if(v > m_timsVmax)
270  {
272  QObject::tr("ERROR in TimsFrame::getVoltageTransformation invalid tims "
273  "calibration, v > m_timsVmax"));
274  }
275  return v;
276 }
277 double
278 TimsFrameBase::getDriftTime(std::size_t scanNum) const
279 {
280  return (m_accumulationTime / (double)m_scanNumber) * ((double)scanNum);
281 }
282 
283 double
285 {
286  return 1 / (m_timsC6 + (m_timsC7 / getVoltageTransformation(scanNum)));
287 }
288 
289 
290 std::size_t
291 TimsFrameBase::getScanNumFromOneOverK0(double one_over_k0) const
292 {
293  double temp = 1 / one_over_k0;
294  temp = temp - m_timsC6;
295  temp = temp / m_timsC7;
296  temp = 1 / temp;
297  temp = temp - m_timsDvStart;
298  temp = temp / m_timsSlope + m_timsTtrans + m_timsNdelay;
299  return (std::size_t)std::round(temp);
300 }
301 
302 bool
304 {
305  if((m_timsDvStart == other.m_timsDvStart) &&
306  (m_timsTtrans == other.m_timsTtrans) &&
307  (m_timsNdelay == other.m_timsNdelay) && (m_timsVmin == other.m_timsVmin) &&
308  (m_timsVmax == other.m_timsVmax) && (m_timsC6 == other.m_timsC6) &&
309  (m_timsC7 == other.m_timsC7) && (m_timsSlope == other.m_timsSlope))
310  {
311  return true;
312  }
313  return false;
314 }
315 
316 
319  std::map<quint32, quint32> &accumulated_scans) const
320 {
321  qDebug();
322  // qDebug();
323  // add flanking peaks
324  pappso::Trace local_trace;
325 
326  MzCalibrationInterface *mz_calibration_p =
328 
329 
330  DataPoint element;
331  for(auto &scan_element : accumulated_scans)
332  {
333  // intensity normalization
334  element.y = ((double)scan_element.second) * 100.0 / m_accumulationTime;
335 
336  // mz calibration
337  element.x = mz_calibration_p->getMzFromTofIndex(scan_element.first);
338 
339  local_trace.push_back(element);
340  }
341  local_trace.sortX();
342 
343  qDebug();
344  // qDebug();
345  return local_trace;
346 }
347 
350  std::map<quint32, quint32> &accumulated_scans) const
351 {
352  qDebug();
353  // qDebug();
354  // add flanking peaks
355  std::vector<quint32> keys;
356  transform(begin(accumulated_scans),
357  end(accumulated_scans),
358  back_inserter(keys),
359  [](std::map<quint32, quint32>::value_type const &pair) {
360  return pair.first;
361  });
362  std::sort(keys.begin(), keys.end());
363  pappso::DataPoint data_point_cumul;
364  data_point_cumul.x = 0;
365  data_point_cumul.y = 0;
366 
367  pappso::Trace local_trace;
368 
369  MzCalibrationInterface *mz_calibration_p =
371 
372 
373  quint32 last_key = 0;
374 
375  for(quint32 key : keys)
376  {
377  if(key == last_key + 1)
378  {
379  // cumulate
380  if(accumulated_scans[key] > accumulated_scans[last_key])
381  {
382  if(data_point_cumul.x == last_key)
383  {
384  // growing peak
385  data_point_cumul.x = key;
386  data_point_cumul.y += accumulated_scans[key];
387  }
388  else
389  {
390  // new peak
391  // flush
392  if(data_point_cumul.y > 0)
393  {
394  // intensity normalization
395  data_point_cumul.y *= 100.0 / m_accumulationTime;
396 
397 
398  // mz calibration
399  data_point_cumul.x =
400  mz_calibration_p->getMzFromTofIndex(data_point_cumul.x);
401  local_trace.push_back(data_point_cumul);
402  }
403 
404  // new point
405  data_point_cumul.x = key;
406  data_point_cumul.y = accumulated_scans[key];
407  }
408  }
409  else
410  {
411  data_point_cumul.y += accumulated_scans[key];
412  }
413  }
414  else
415  {
416  // flush
417  if(data_point_cumul.y > 0)
418  {
419  // intensity normalization
420  data_point_cumul.y *= 100.0 / m_accumulationTime;
421 
422 
423  // qDebug() << "raw data x=" << data_point_cumul.x;
424  // mz calibration
425  data_point_cumul.x =
426  mz_calibration_p->getMzFromTofIndex(data_point_cumul.x);
427  // qDebug() << "mz=" << data_point_cumul.x;
428  local_trace.push_back(data_point_cumul);
429  }
430 
431  // new point
432  data_point_cumul.x = key;
433  data_point_cumul.y = accumulated_scans[key];
434  }
435 
436  last_key = key;
437  }
438  // flush
439  if(data_point_cumul.y > 0)
440  {
441  // intensity normalization
442  data_point_cumul.y *= 100.0 / m_accumulationTime;
443 
444 
445  // mz calibration
446  data_point_cumul.x =
447  mz_calibration_p->getMzFromTofIndex(data_point_cumul.x);
448  local_trace.push_back(data_point_cumul);
449  }
450 
451  local_trace.sortX();
452  qDebug();
453  // qDebug();
454  return local_trace;
455 }
456 
459 {
460  if(msp_mzCalibration == nullptr)
461  {
462 
464  QObject::tr("ERROR in %1, %2, %3 msp_mzCalibration is null")
465  .arg(__FILE__)
466  .arg(__FUNCTION__)
467  .arg(__LINE__));
468  }
469  return msp_mzCalibration;
470 }
471 
472 void
474  MzCalibrationInterfaceSPtr mzCalibration)
475 {
476 
477  if(mzCalibration == nullptr)
478  {
479 
481  QObject::tr("ERROR in %1, %2, %3 msp_mzCalibration is null")
482  .arg(__FILE__)
483  .arg(__FUNCTION__)
484  .arg(__LINE__));
485  }
486  msp_mzCalibration = mzCalibration;
487 }
488 
489 
490 quint32
492 {
493  quint32 max_value = 0;
494  for(quint32 i = 0; i < m_scanNumber; i++)
495  {
496  qDebug() << "m_scanNumber=" << m_scanNumber << " i=" << i;
497  std::vector<quint32> index_list = getScanIndexList(i);
498  auto it = std::max_element(index_list.begin(), index_list.end());
499  if(it != index_list.end())
500  {
501  max_value = std::max(max_value, *it);
502  }
503  }
504  return max_value;
505 }
506 
507 std::vector<quint32>
508 TimsFrameBase::getScanIndexList(std::size_t scanNum) const
509 {
510  throw PappsoException(
511  QObject::tr(
512  "ERROR unable to getScanIndexList in TimsFrameBase for scan number %1")
513  .arg(scanNum));
514 }
515 
516 
517 std::vector<quint32>
518 TimsFrameBase::getScanIntensities(std::size_t scanNum) const
519 {
520  throw PappsoException(
521  QObject::tr(
522  "ERROR unable to getScanIntensities in TimsFrameBase for scan number %1")
523  .arg(scanNum));
524 }
525 
526 Trace
528  std::size_t mz_index_lower_bound,
529  std::size_t mz_index_upper_bound,
530  XicExtractMethod method) const
531 {
532  Trace im_trace;
533  DataPoint data_point;
534  for(quint32 i = 0; i < m_scanNumber; i++)
535  {
536  data_point.x = i;
537  data_point.y = 0;
538  qDebug() << "m_scanNumber=" << m_scanNumber << " i=" << i;
539  std::vector<quint32> index_list = getScanIndexList(i);
540  auto it_lower = std::find_if(index_list.begin(),
541  index_list.end(),
542  [mz_index_lower_bound](quint32 to_compare) {
543  if(to_compare < mz_index_lower_bound)
544  {
545  return false;
546  }
547  return true;
548  });
549 
550 
551  if(it_lower == index_list.end())
552  {
553  }
554  else
555  {
556 
557 
558  auto it_upper =
559  std::find_if(index_list.begin(),
560  index_list.end(),
561  [mz_index_upper_bound](quint32 to_compare) {
562  if(mz_index_upper_bound >= to_compare)
563  {
564  return false;
565  }
566  return true;
567  });
568  std::vector<quint32> intensity_list = getScanIntensities(i);
569  for(int j = std::distance(index_list.begin(), it_lower);
570  j < std::distance(index_list.begin(), it_upper);
571  j++)
572  {
573  if(method == XicExtractMethod::sum)
574  {
575  data_point.y += intensity_list[j];
576  }
577  else
578  {
579  data_point.y =
580  std::max((double)intensity_list[j], data_point.y);
581  }
582  }
583  }
584  im_trace.push_back(data_point);
585  }
586  qDebug();
587  return im_trace;
588 }
589 // namespace pappso
590 } // namespace pappso
virtual double getMzFromTofIndex(quint32 tof_index)=0
get m/z from time of flight raw index
double m_accumulationTime
accumulation time in milliseconds
double getTime() const
virtual quint64 cumulateSingleScanIntensities(std::size_t scanNum) const
double getVoltageTransformation(std::size_t scanNum) const
get voltage for a given scan number
virtual Trace cumulateScanToTrace(std::size_t scanNumBegin, std::size_t scanNumEnd) const
cumulate spectrum given a scan number range need the binary file The intensities are normalized with ...
virtual std::size_t getTotalNumberOfScans() const
get the number of scans contained in this frame each scan represents an ion mobility slice
virtual Trace getIonMobilityTraceByMzIndexRange(std::size_t mz_index_lower_bound, std::size_t mz_index_upper_bound, XicExtractMethod method) const
get a mobility trace cumulating intensities inside the given mass index range
virtual std::size_t getNbrPeaks(std::size_t scanNum) const
get the number of peaks in this spectrum need the binary file
MzCalibrationInterfaceSPtr msp_mzCalibration
virtual MassSpectrumSPtr getMassSpectrumSPtr(std::size_t scanNum) const
get Mass spectrum with peaks for this scan number need the binary file
TimsFrameBase(std::size_t timsId, quint32 scanNum)
constructor for binary independant tims frame
virtual std::vector< quint32 > getScanIndexList(std::size_t scanNum) const
get raw index list for one given scan index are not TOF nor m/z, just index on digitizer
double getDriftTime(std::size_t scanNum) const
get drift time of a scan number in milliseconds
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
double m_time
retention time
virtual quint64 cumulateScansIntensities(std::size_t scanNumBegin, std::size_t scanNumEnd) const
void setAccumulationTime(double accumulation_time_ms)
quint32 m_scanNumber
total number of scans contained in this frame
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 The intensities are NOT normalized with respe...
void setTime(double time)
std::size_t m_timsId
Tims frame database id (the SQL identifier of this frame)
pappso::Trace getTraceFromCumulatedScans(std::map< quint32, quint32 > &accumulated_scans) const
transform accumulation of raw scans into a real mass spectrum
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...
virtual quint32 getMaximumRawMassIndex() const
get the maximum raw mass index contained in this frame
unsigned int getMsLevel() const
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)
virtual const MzCalibrationInterfaceSPtr & getMzCalibrationInterfaceSPtr() const final
get the MzCalibration model to compute mz and TOF for this frame
std::size_t getScanNumFromOneOverK0(double one_over_k0) const
get the scan number from a given 1/Ko mobility value
void setMsMsType(quint8 type)
void setMzCalibration(double T1_frame, double T2_frame, double digitizerTimebase, double digitizerDelay, double C0, double C1, double C2, double C3, double C4, double T1_ref, double T2_ref, double dC1, double dC2)
double getOneOverK0Transformation(std::size_t scanNum) const
get 1/K0 value of a given scan (mobility value)
bool checkScanNum(std::size_t scanNum) const
check that this scan number exists
void setMzCalibrationInterfaceSPtr(MzCalibrationInterfaceSPtr mzCalibration)
std::size_t getId() const
virtual std::vector< quint32 > getScanIntensities(std::size_t scanNum) const
get raw intensities without transformation from one scan it needs intensity normalization
A simple container of DataPoint instances.
Definition: trace.h:148
void sortX()
Definition: trace.cpp:956
implement Bruker's model type 1 formula to compute m/z
tries to keep as much as possible monoisotopes, removing any possible C13 peaks and changes multichar...
Definition: aa.cpp:39
std::shared_ptr< MzCalibrationInterface > MzCalibrationInterfaceSPtr
std::shared_ptr< MassSpectrum > MassSpectrumSPtr
Definition: massspectrum.h:54
XicExtractMethod
Definition: types.h:201
@ sum
sum of intensities
pappso_double x
Definition: datapoint.h:23
pappso_double y
Definition: datapoint.h:24
handle a single Bruker's TimsTof frame without binary data