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
35namespace pappso
36{
37
38TimsFrameBase::TimsFrameBase(std::size_t timsId, quint32 scanNum)
39{
40 qDebug() << timsId;
41 m_timsId = timsId;
42
43 m_scanNumber = scanNum;
44}
45
46TimsFrameBase::TimsFrameBase([[maybe_unused]] const TimsFrameBase &other)
47{
48}
49
51{
52}
53
54void
55TimsFrameBase::setAccumulationTime(double accumulation_time_ms)
56{
57 m_accumulationTime = accumulation_time_ms;
58}
59
60
61void
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
101bool
102TimsFrameBase::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
115std::size_t
116TimsFrameBase::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
124std::size_t
126{
127 return m_scanNumber;
128}
129
131TimsFrameBase::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}
138Trace
139TimsFrameBase::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}
148void
149TimsFrameBase::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}
161void
163{
164 m_time = time;
165}
166
167void
169{
170
171 qDebug() << " m_msMsType=" << type;
172 m_msMsType = type;
173}
174
175unsigned int
177{
178 if(m_msMsType == 0)
179 return 1;
180 return 2;
181}
182
183double
185{
186 return m_time;
187}
188
189std::size_t
191{
192 return m_timsId;
193}
194void
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
222 (C3 - m_timsDvStart) / C1; // //C3 from TimsCalibration // C2 from
223 // TimsCalibration // C1 from TimsCalibration
224}
225double
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}
247double
248TimsFrameBase::getDriftTime(std::size_t scanNum) const
249{
250 return (m_accumulationTime / (double)m_scanNumber) * ((double)scanNum);
251}
252
253double
255{
256 return 1 / (m_timsC6 + (m_timsC7 / getVoltageTransformation(scanNum)));
257}
258
259
260std::size_t
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
272bool
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
442void
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
460quint32
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
477std::vector<quint32>
478TimsFrameBase::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
487std::vector<quint32>
488TimsFrameBase::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
496Trace
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