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 void
150 TimsFrameBase::cumulateScansInRawMap(std::map<quint32, quint32> &rawSpectrum
151  [[maybe_unused]],
152  std::size_t scanNumBegin,
153  std::size_t scanNumEnd) const
154 {
155  throw PappsoException(
156  QObject::tr(
157  "ERROR unable to cumulateScansInRawMap in TimsFrameBase for scan "
158  "number begin %1 end %2")
159  .arg(scanNumBegin)
160  .arg(scanNumEnd));
161 }
162 
163 
164 quint64
166 {
167  throw PappsoException(
168  QObject::tr(
169  "ERROR unable to cumulateSingleScanIntensities in TimsFrameBase for scan "
170  "number %1.").arg(scanNum));
171 
172  return 0;
173 }
174 
175 
176 quint64
178  std::size_t scanNumEnd) const
179 {
180  throw PappsoException(
181  QObject::tr(
182  "ERROR unable to cumulateScansInRawMap in TimsFrameBase for scan "
183  "number begin %1 end %2")
184  .arg(scanNumBegin)
185  .arg(scanNumEnd));
186 
187  return 0;
188 }
189 
190 void
192 {
193  m_time = time;
194 }
195 
196 void
198 {
199 
200  qDebug() << " m_msMsType=" << type;
201  m_msMsType = type;
202 }
203 
204 unsigned int
206 {
207  if(m_msMsType == 0)
208  return 1;
209  return 2;
210 }
211 
212 double
214 {
215  return m_time;
216 }
217 
218 std::size_t
220 {
221  return m_timsId;
222 }
223 void
225  double C0,
226  double C1,
227  double C2,
228  double C3,
229  double C4,
230  [[maybe_unused]] double C5,
231  double C6,
232  double C7,
233  double C8,
234  double C9)
235 {
236  if(tims_model_type != 2)
237  {
238  throw pappso::PappsoException(QObject::tr(
239  "ERROR in TimsFrame::setTimsCalibration tims_model_type != 2"));
240  }
241  m_timsDvStart = C2; // C2 from TimsCalibration
242  m_timsTtrans = C4; // C4 from TimsCalibration
243  m_timsNdelay = C0; // C0 from TimsCalibration
244  m_timsVmin = C8; // C8 from TimsCalibration
245  m_timsVmax = C9; // C9 from TimsCalibration
246  m_timsC6 = C6;
247  m_timsC7 = C7;
248 
249 
250  m_timsSlope =
251  (C3 - m_timsDvStart) / C1; // //C3 from TimsCalibration // C2 from
252  // TimsCalibration // C1 from TimsCalibration
253 }
254 double
255 TimsFrameBase::getVoltageTransformation(std::size_t scanNum) const
256 {
257  double v = m_timsDvStart +
258  m_timsSlope * ((double)scanNum - m_timsTtrans - m_timsNdelay);
259 
260  if(v < m_timsVmin)
261  {
263  QObject::tr("ERROR in TimsFrame::getVoltageTransformation invalid tims "
264  "calibration, v < m_timsVmin"));
265  }
266 
267 
268  if(v > m_timsVmax)
269  {
271  QObject::tr("ERROR in TimsFrame::getVoltageTransformation invalid tims "
272  "calibration, v > m_timsVmax"));
273  }
274  return v;
275 }
276 double
277 TimsFrameBase::getDriftTime(std::size_t scanNum) const
278 {
279  return (m_accumulationTime / (double)m_scanNumber) * ((double)scanNum);
280 }
281 
282 double
284 {
285  return 1 / (m_timsC6 + (m_timsC7 / getVoltageTransformation(scanNum)));
286 }
287 
288 
289 std::size_t
290 TimsFrameBase::getScanNumFromOneOverK0(double one_over_k0) const
291 {
292  double temp = 1 / one_over_k0;
293  temp = temp - m_timsC6;
294  temp = temp / m_timsC7;
295  temp = 1 / temp;
296  temp = temp - m_timsDvStart;
297  temp = temp / m_timsSlope + m_timsTtrans + m_timsNdelay;
298  return (std::size_t)std::round(temp);
299 }
300 
301 bool
303 {
304  if((m_timsDvStart == other.m_timsDvStart) &&
305  (m_timsTtrans == other.m_timsTtrans) &&
306  (m_timsNdelay == other.m_timsNdelay) && (m_timsVmin == other.m_timsVmin) &&
307  (m_timsVmax == other.m_timsVmax) && (m_timsC6 == other.m_timsC6) &&
308  (m_timsC7 == other.m_timsC7) && (m_timsSlope == other.m_timsSlope))
309  {
310  return true;
311  }
312  return false;
313 }
314 
315 
318  std::map<quint32, quint32> &accumulated_scans) const
319 {
320  qDebug();
321  // qDebug();
322  // add flanking peaks
323  pappso::Trace local_trace;
324 
325  MzCalibrationInterface *mz_calibration_p =
327 
328 
329  DataPoint element;
330  for(auto &scan_element : accumulated_scans)
331  {
332  // intensity normalization
333  element.y = ((double)scan_element.second) * 100.0 / m_accumulationTime;
334 
335  // mz calibration
336  element.x = mz_calibration_p->getMzFromTofIndex(scan_element.first);
337 
338  local_trace.push_back(element);
339  }
340  local_trace.sortX();
341 
342  qDebug();
343  // qDebug();
344  return local_trace;
345 }
346 
349  std::map<quint32, quint32> &accumulated_scans) const
350 {
351  qDebug();
352  // qDebug();
353  // add flanking peaks
354  std::vector<quint32> keys;
355  transform(begin(accumulated_scans),
356  end(accumulated_scans),
357  back_inserter(keys),
358  [](std::map<quint32, quint32>::value_type const &pair) {
359  return pair.first;
360  });
361  std::sort(keys.begin(), keys.end());
362  pappso::DataPoint data_point_cumul;
363  data_point_cumul.x = 0;
364  data_point_cumul.y = 0;
365 
366  pappso::Trace local_trace;
367 
368  MzCalibrationInterface *mz_calibration_p =
370 
371 
372  quint32 last_key = 0;
373 
374  for(quint32 key : keys)
375  {
376  if(key == last_key + 1)
377  {
378  // cumulate
379  if(accumulated_scans[key] > accumulated_scans[last_key])
380  {
381  if(data_point_cumul.x == last_key)
382  {
383  // growing peak
384  data_point_cumul.x = key;
385  data_point_cumul.y += accumulated_scans[key];
386  }
387  else
388  {
389  // new peak
390  // flush
391  if(data_point_cumul.y > 0)
392  {
393  // intensity normalization
394  data_point_cumul.y *= 100.0 / m_accumulationTime;
395 
396 
397  // mz calibration
398  data_point_cumul.x =
399  mz_calibration_p->getMzFromTofIndex(data_point_cumul.x);
400  local_trace.push_back(data_point_cumul);
401  }
402 
403  // new point
404  data_point_cumul.x = key;
405  data_point_cumul.y = accumulated_scans[key];
406  }
407  }
408  else
409  {
410  data_point_cumul.y += accumulated_scans[key];
411  }
412  }
413  else
414  {
415  // flush
416  if(data_point_cumul.y > 0)
417  {
418  // intensity normalization
419  data_point_cumul.y *= 100.0 / m_accumulationTime;
420 
421 
422  // qDebug() << "raw data x=" << data_point_cumul.x;
423  // mz calibration
424  data_point_cumul.x =
425  mz_calibration_p->getMzFromTofIndex(data_point_cumul.x);
426  // qDebug() << "mz=" << data_point_cumul.x;
427  local_trace.push_back(data_point_cumul);
428  }
429 
430  // new point
431  data_point_cumul.x = key;
432  data_point_cumul.y = accumulated_scans[key];
433  }
434 
435  last_key = key;
436  }
437  // flush
438  if(data_point_cumul.y > 0)
439  {
440  // intensity normalization
441  data_point_cumul.y *= 100.0 / m_accumulationTime;
442 
443 
444  // mz calibration
445  data_point_cumul.x =
446  mz_calibration_p->getMzFromTofIndex(data_point_cumul.x);
447  local_trace.push_back(data_point_cumul);
448  }
449 
450  local_trace.sortX();
451  qDebug();
452  // qDebug();
453  return local_trace;
454 }
455 
458 {
459  if(msp_mzCalibration == nullptr)
460  {
461 
463  QObject::tr("ERROR in %1, %2, %3 msp_mzCalibration is null")
464  .arg(__FILE__)
465  .arg(__FUNCTION__)
466  .arg(__LINE__));
467  }
468  return msp_mzCalibration;
469 }
470 
471 void
473  MzCalibrationInterfaceSPtr mzCalibration)
474 {
475 
476  if(mzCalibration == nullptr)
477  {
478 
480  QObject::tr("ERROR in %1, %2, %3 msp_mzCalibration is null")
481  .arg(__FILE__)
482  .arg(__FUNCTION__)
483  .arg(__LINE__));
484  }
485  msp_mzCalibration = mzCalibration;
486 }
487 
488 
489 quint32
491 {
492  quint32 max_value = 0;
493  for(quint32 i = 0; i < m_scanNumber; i++)
494  {
495  qDebug() << "m_scanNumber=" << m_scanNumber << " i=" << i;
496  std::vector<quint32> index_list = getScanIndexList(i);
497  auto it = std::max_element(index_list.begin(), index_list.end());
498  if(it != index_list.end())
499  {
500  max_value = std::max(max_value, *it);
501  }
502  }
503  return max_value;
504 }
505 
506 std::vector<quint32>
507 TimsFrameBase::getScanIndexList(std::size_t scanNum) const
508 {
509  throw PappsoException(
510  QObject::tr(
511  "ERROR unable to getScanIndexList in TimsFrameBase for scan number %1")
512  .arg(scanNum));
513 }
514 
515 
516 std::vector<quint32>
517 TimsFrameBase::getScanIntensities(std::size_t scanNum) const
518 {
519  throw PappsoException(
520  QObject::tr(
521  "ERROR unable to getScanIntensities in TimsFrameBase for scan number %1")
522  .arg(scanNum));
523 }
524 
525 Trace
527  std::size_t mz_index_lower_bound,
528  std::size_t mz_index_upper_bound,
529  XicExtractMethod method) const
530 {
531  Trace im_trace;
532  DataPoint data_point;
533  for(quint32 i = 0; i < m_scanNumber; i++)
534  {
535  data_point.x = i;
536  data_point.y = 0;
537  qDebug() << "m_scanNumber=" << m_scanNumber << " i=" << i;
538  std::vector<quint32> index_list = getScanIndexList(i);
539  auto it_lower = std::find_if(index_list.begin(),
540  index_list.end(),
541  [mz_index_lower_bound](quint32 to_compare) {
542  if(to_compare < mz_index_lower_bound)
543  {
544  return false;
545  }
546  return true;
547  });
548 
549 
550  if(it_lower == index_list.end())
551  {
552  }
553  else
554  {
555 
556 
557  auto it_upper =
558  std::find_if(index_list.begin(),
559  index_list.end(),
560  [mz_index_upper_bound](quint32 to_compare) {
561  if(mz_index_upper_bound >= to_compare)
562  {
563  return false;
564  }
565  return true;
566  });
567  std::vector<quint32> intensity_list = getScanIntensities(i);
568  for(int j = std::distance(index_list.begin(), it_lower);
569  j < std::distance(index_list.begin(), it_upper);
570  j++)
571  {
572  if(method == XicExtractMethod::sum)
573  {
574  data_point.y += intensity_list[j];
575  }
576  else
577  {
578  data_point.y =
579  std::max((double)intensity_list[j], data_point.y);
580  }
581  }
582  }
583  im_trace.push_back(data_point);
584  }
585  qDebug();
586  return im_trace;
587 }
588 // namespace pappso
589 } // 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:38
Q_INVOKABLE void sortX()
Definition: trace.cpp:983
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