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 <cmath>
33 #include <algorithm>
34 
35 namespace pappso
36 {
37 
38 TimsFrameBase::TimsFrameBase(std::size_t timsId, quint32 scanNum)
39 {
40  qDebug() << timsId;
41  m_timsId = timsId;
42 
43  m_scanNumber = scanNum;
44 }
45 
46 TimsFrameBase::TimsFrameBase([[maybe_unused]] const TimsFrameBase &other)
47 {
48 }
49 
51 {
52 }
53 
54 void
55 TimsFrameBase::setAccumulationTime(double accumulation_time_ms)
56 {
57  m_accumulationTime = accumulation_time_ms;
58 }
59 
60 
61 void
63  double T2_frame,
64  double digitizerTimebase,
65  double digitizerDelay,
66  double C0,
67  double C1,
68  double C2,
69  double C3,
70  double C4,
71  double T1_ref,
72  double T2_ref,
73  double dC1,
74  double dC2)
75 {
76 
77  /* MzCalibrationModel1 mzCalibration(temperature_correction,
78  digitizerTimebase,
79  digitizerDelay,
80  C0,
81  C1,
82  C2,
83  C3,
84  C4);
85  */
86  msp_mzCalibration = std::make_shared<MzCalibrationModel1>(T1_frame,
87  T2_frame,
88  digitizerTimebase,
89  digitizerDelay,
90  C0,
91  C1,
92  C2,
93  C3,
94  C4,
95  T1_ref,
96  T2_ref,
97  dC1,
98  dC2);
99 }
100 
101 bool
102 TimsFrameBase::checkScanNum(std::size_t scanNum) const
103 {
104  if(scanNum >= m_scanNumber)
105  {
107  QObject::tr("Invalid scan number : scanNum %1 > m_scanNumber %2")
108  .arg(scanNum)
109  .arg(m_scanNumber));
110  }
111 
112  return true;
113 }
114 
115 std::size_t
116 TimsFrameBase::getNbrPeaks(std::size_t scanNum) const
117 {
118  throw PappsoException(
119  QObject::tr(
120  "ERROR unable to get number of peaks in TimsFrameBase for scan number %1")
121  .arg(scanNum));
122 }
123 
124 std::size_t
126 {
127  return m_scanNumber;
128 }
129 
131 TimsFrameBase::getMassSpectrumSPtr(std::size_t scanNum) const
132 {
133  throw PappsoException(
134  QObject::tr(
135  "ERROR unable to getMassSpectrumSPtr in TimsFrameBase for scan number %1")
136  .arg(scanNum));
137 }
138 Trace
139 TimsFrameBase::cumulateScanToTrace(std::size_t scanNumBegin,
140  std::size_t scanNumEnd) const
141 {
142  throw PappsoException(
143  QObject::tr("ERROR unable to cumulateScanToTrace in TimsFrameBase for scan "
144  "number begin %1 end %2")
145  .arg(scanNumBegin)
146  .arg(scanNumEnd));
147 }
148 void
149 TimsFrameBase::cumulateScansInRawMap(std::map<quint32, quint32> &rawSpectrum
150  [[maybe_unused]],
151  std::size_t scanNumBegin,
152  std::size_t scanNumEnd) const
153 {
154  throw PappsoException(
155  QObject::tr(
156  "ERROR unable to cumulateScansInRawMap in TimsFrameBase for scan "
157  "number begin %1 end %2")
158  .arg(scanNumBegin)
159  .arg(scanNumEnd));
160 }
161 void
163 {
164  m_time = time;
165 }
166 
167 void
169 {
170 
171  qDebug() << " m_msMsType=" << type;
172  m_msMsType = type;
173 }
174 
175 unsigned int
177 {
178  if(m_msMsType == 0)
179  return 1;
180  return 2;
181 }
182 
183 double
185 {
186  return m_time;
187 }
188 
189 std::size_t
191 {
192  return m_timsId;
193 }
194 void
196  double C0,
197  double C1,
198  double C2,
199  double C3,
200  double C4,
201  [[maybe_unused]] double C5,
202  double C6,
203  double C7,
204  double C8,
205  double C9)
206 {
207  if(tims_model_type != 2)
208  {
209  throw pappso::PappsoException(QObject::tr(
210  "ERROR in TimsFrame::setTimsCalibration tims_model_type != 2"));
211  }
212  m_timsDvStart = C2; // C2 from TimsCalibration
213  m_timsTtrans = C4; // C4 from TimsCalibration
214  m_timsNdelay = C0; // C0 from TimsCalibration
215  m_timsVmin = C8; // C8 from TimsCalibration
216  m_timsVmax = C9; // C9 from TimsCalibration
217  m_timsC6 = C6;
218  m_timsC7 = C7;
219 
220 
221  m_timsSlope =
222  (C3 - m_timsDvStart) / C1; // //C3 from TimsCalibration // C2 from
223  // TimsCalibration // C1 from TimsCalibration
224 }
225 double
226 TimsFrameBase::getVoltageTransformation(std::size_t scanNum) const
227 {
228  double v = m_timsDvStart +
229  m_timsSlope * ((double)scanNum - m_timsTtrans - m_timsNdelay);
230 
231  if(v < m_timsVmin)
232  {
234  QObject::tr("ERROR in TimsFrame::getVoltageTransformation invalid tims "
235  "calibration, v < m_timsVmin"));
236  }
237 
238 
239  if(v > m_timsVmax)
240  {
242  QObject::tr("ERROR in TimsFrame::getVoltageTransformation invalid tims "
243  "calibration, v > m_timsVmax"));
244  }
245  return v;
246 }
247 double
248 TimsFrameBase::getDriftTime(std::size_t scanNum) const
249 {
250  return (m_accumulationTime / (double)m_scanNumber) * ((double)scanNum);
251 }
252 
253 double
255 {
256  return 1 / (m_timsC6 + (m_timsC7 / getVoltageTransformation(scanNum)));
257 }
258 
259 
260 std::size_t
261 TimsFrameBase::getScanNumFromOneOverK0(double one_over_k0) const
262 {
263  double temp = 1 / one_over_k0;
264  temp = temp - m_timsC6;
265  temp = temp / m_timsC7;
266  temp = 1 / temp;
267  temp = temp - m_timsDvStart;
268  temp = temp / m_timsSlope + m_timsTtrans + m_timsNdelay;
269  return (std::size_t)std::round(temp);
270 }
271 
272 bool
274 {
275  if((m_timsDvStart == other.m_timsDvStart) &&
276  (m_timsTtrans == other.m_timsTtrans) &&
277  (m_timsNdelay == other.m_timsNdelay) && (m_timsVmin == other.m_timsVmin) &&
278  (m_timsVmax == other.m_timsVmax) && (m_timsC6 == other.m_timsC6) &&
279  (m_timsC7 == other.m_timsC7) && (m_timsSlope == other.m_timsSlope))
280  {
281  return true;
282  }
283  return false;
284 }
285 
286 
289  std::map<quint32, quint32> &accumulated_scans) const
290 {
291  qDebug();
292  // qDebug();
293  // add flanking peaks
294  pappso::Trace local_trace;
295 
296  MzCalibrationInterface *mz_calibration_p =
298 
299 
300  DataPoint element;
301  for(auto &scan_element : accumulated_scans)
302  {
303  // intensity normalization
304  element.y = ((double)scan_element.second) * 100.0 / m_accumulationTime;
305 
306  // mz calibration
307  element.x = mz_calibration_p->getMzFromTofIndex(scan_element.first);
308 
309  local_trace.push_back(element);
310  }
311  local_trace.sortX();
312 
313  qDebug();
314  // qDebug();
315  return local_trace;
316 }
317 
320  std::map<quint32, quint32> &accumulated_scans) const
321 {
322  qDebug();
323  // qDebug();
324  // add flanking peaks
325  std::vector<quint32> keys;
326  transform(begin(accumulated_scans),
327  end(accumulated_scans),
328  back_inserter(keys),
329  [](std::map<quint32, quint32>::value_type const &pair) {
330  return pair.first;
331  });
332  std::sort(keys.begin(), keys.end());
333  pappso::DataPoint data_point_cumul;
334  data_point_cumul.x = 0;
335  data_point_cumul.y = 0;
336 
337  pappso::Trace local_trace;
338 
339  MzCalibrationInterface *mz_calibration_p =
341 
342 
343  quint32 last_key = 0;
344 
345  for(quint32 key : keys)
346  {
347  if(key == last_key + 1)
348  {
349  // cumulate
350  if(accumulated_scans[key] > accumulated_scans[last_key])
351  {
352  if(data_point_cumul.x == last_key)
353  {
354  // growing peak
355  data_point_cumul.x = key;
356  data_point_cumul.y += accumulated_scans[key];
357  }
358  else
359  {
360  // new peak
361  // flush
362  if(data_point_cumul.y > 0)
363  {
364  // intensity normalization
365  data_point_cumul.y *= 100.0 / m_accumulationTime;
366 
367 
368  // mz calibration
369  data_point_cumul.x =
370  mz_calibration_p->getMzFromTofIndex(data_point_cumul.x);
371  local_trace.push_back(data_point_cumul);
372  }
373 
374  // new point
375  data_point_cumul.x = key;
376  data_point_cumul.y = accumulated_scans[key];
377  }
378  }
379  else
380  {
381  data_point_cumul.y += accumulated_scans[key];
382  }
383  }
384  else
385  {
386  // flush
387  if(data_point_cumul.y > 0)
388  {
389  // intensity normalization
390  data_point_cumul.y *= 100.0 / m_accumulationTime;
391 
392 
393  // qDebug() << "raw data x=" << data_point_cumul.x;
394  // mz calibration
395  data_point_cumul.x =
396  mz_calibration_p->getMzFromTofIndex(data_point_cumul.x);
397  // qDebug() << "mz=" << data_point_cumul.x;
398  local_trace.push_back(data_point_cumul);
399  }
400 
401  // new point
402  data_point_cumul.x = key;
403  data_point_cumul.y = accumulated_scans[key];
404  }
405 
406  last_key = key;
407  }
408  // flush
409  if(data_point_cumul.y > 0)
410  {
411  // intensity normalization
412  data_point_cumul.y *= 100.0 / m_accumulationTime;
413 
414 
415  // mz calibration
416  data_point_cumul.x =
417  mz_calibration_p->getMzFromTofIndex(data_point_cumul.x);
418  local_trace.push_back(data_point_cumul);
419  }
420 
421  local_trace.sortX();
422  qDebug();
423  // qDebug();
424  return local_trace;
425 }
426 
429 {
430  if(msp_mzCalibration == nullptr)
431  {
432 
434  QObject::tr("ERROR in %1, %2, %3 msp_mzCalibration is null")
435  .arg(__FILE__)
436  .arg(__FUNCTION__)
437  .arg(__LINE__));
438  }
439  return msp_mzCalibration;
440 }
441 
442 void
444  MzCalibrationInterfaceSPtr mzCalibration)
445 {
446 
447  if(mzCalibration == nullptr)
448  {
449 
451  QObject::tr("ERROR in %1, %2, %3 msp_mzCalibration is null")
452  .arg(__FILE__)
453  .arg(__FUNCTION__)
454  .arg(__LINE__));
455  }
456  msp_mzCalibration = mzCalibration;
457 }
458 
459 
460 quint32
462 {
463  quint32 max_value = 0;
464  for(quint32 i = 0; i < m_scanNumber; i++)
465  {
466  qDebug() << "m_scanNumber=" << m_scanNumber << " i=" << i;
467  std::vector<quint32> index_list = getScanIndexList(i);
468  auto it = std::max_element(index_list.begin(), index_list.end());
469  if(it != index_list.end())
470  {
471  max_value = std::max(max_value, *it);
472  }
473  }
474  return max_value;
475 }
476 
477 std::vector<quint32>
478 TimsFrameBase::getScanIndexList(std::size_t scanNum) const
479 {
480  throw PappsoException(
481  QObject::tr(
482  "ERROR unable to getScanIndexList in TimsFrameBase for scan number %1")
483  .arg(scanNum));
484 }
485 
486 
487 std::vector<quint32>
488 TimsFrameBase::getScanIntensities(std::size_t scanNum) const
489 {
490  throw PappsoException(
491  QObject::tr(
492  "ERROR unable to getScanIntensities in TimsFrameBase for scan number %1")
493  .arg(scanNum));
494 }
495 
496 Trace
498  std::size_t mz_index_lower_bound,
499  std::size_t mz_index_upper_bound,
500  XicExtractMethod method) const
501 {
502  Trace im_trace;
503  DataPoint data_point;
504  for(quint32 i = 0; i < m_scanNumber; i++)
505  {
506  data_point.x = i;
507  data_point.y = 0;
508  qDebug() << "m_scanNumber=" << m_scanNumber << " i=" << i;
509  std::vector<quint32> index_list = getScanIndexList(i);
510  auto it_lower = std::find_if(index_list.begin(),
511  index_list.end(),
512  [mz_index_lower_bound](quint32 to_compare) {
513  if(to_compare < mz_index_lower_bound)
514  {
515  return false;
516  }
517  return true;
518  });
519 
520 
521  if(it_lower == index_list.end())
522  {
523  }
524  else
525  {
526 
527 
528  auto it_upper =
529  std::find_if(index_list.begin(),
530  index_list.end(),
531  [mz_index_upper_bound](quint32 to_compare) {
532  if(mz_index_upper_bound >= to_compare)
533  {
534  return false;
535  }
536  return true;
537  });
538  std::vector<quint32> intensity_list = getScanIntensities(i);
539  for(int j = std::distance(index_list.begin(), it_lower);
540  j < std::distance(index_list.begin(), it_upper);
541  j++)
542  {
543  if(method == XicExtractMethod::sum)
544  {
545  data_point.y += intensity_list[j];
546  }
547  else
548  {
549  data_point.y =
550  std::max((double)intensity_list[j], data_point.y);
551  }
552  }
553  }
554  im_trace.push_back(data_point);
555  }
556  qDebug();
557  return im_trace;
558 }
559 // namespace pappso
560 } // 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
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
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
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
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:147
void sortX()
Definition: trace.cpp:936
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:200
@ sum
sum of intensities
pappso_double x
Definition: datapoint.h:22
pappso_double y
Definition: datapoint.h:23
handle a single Bruker's TimsTof frame without binary data