libpappsomspp
Library for mass spectrometry
Loading...
Searching...
No Matches
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
36namespace pappso
37{
38
39TimsFrameBase::TimsFrameBase(std::size_t timsId, quint32 scanNum)
40{
41 qDebug() << timsId;
42 m_timsId = timsId;
43
44 m_scanNumber = scanNum;
45}
46
47TimsFrameBase::TimsFrameBase([[maybe_unused]] const TimsFrameBase &other)
48{
49}
50
54
55void
56TimsFrameBase::setAccumulationTime(double accumulation_time_ms)
57{
58 m_accumulationTime = accumulation_time_ms;
59}
60
61
62void
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
102bool
103TimsFrameBase::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
116std::size_t
117TimsFrameBase::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
125std::size_t
130
132TimsFrameBase::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}
139Trace
140TimsFrameBase::cumulateScansToTrace(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
150Trace
152 std::size_t mzindex_merge_window [[maybe_unused]],
153 std::size_t scanNumBegin [[maybe_unused]],
154 std::size_t scanNumEnd [[maybe_unused]],
155 quint32 &minimum_index [[maybe_unused]],
156 quint32 &maximum_index [[maybe_unused]]) const
157{
158 throw PappsoException(QObject::tr("Non implemented function %1 %2 %3")
159 .arg(__FILE__)
160 .arg(__FUNCTION__)
161 .arg(__LINE__));
162}
163
164void
165TimsFrameBase::cumulateScansInRawMap(std::map<quint32, quint32> &rawSpectrum
166 [[maybe_unused]],
167 std::size_t scanNumBegin,
168 std::size_t scanNumEnd) const
169{
170 throw PappsoException(
171 QObject::tr(
172 "ERROR unable to cumulateScansInRawMap in TimsFrameBase for scan "
173 "number begin %1 end %2")
174 .arg(scanNumBegin)
175 .arg(scanNumEnd));
176}
177
178
179quint64
181{
182 throw PappsoException(
183 QObject::tr(
184 "ERROR unable to cumulateSingleScanIntensities in TimsFrameBase for scan "
185 "number %1.")
186 .arg(scanNum));
187
188 return 0;
189}
190
191
192quint64
194 std::size_t scanNumEnd) const
195{
196 throw PappsoException(
197 QObject::tr(
198 "ERROR unable to cumulateScansInRawMap in TimsFrameBase for scan "
199 "number begin %1 end %2")
200 .arg(scanNumBegin)
201 .arg(scanNumEnd));
202
203 return 0;
204}
205
206void
208{
209 m_time = time;
210}
211
212void
214{
215
216 qDebug() << " m_msMsType=" << type;
217 m_msMsType = type;
218}
219
220unsigned int
222{
223 if(m_msMsType == 0)
224 return 1;
225 return 2;
226}
227
228double
230{
231 return m_time;
232}
233
234std::size_t
236{
237 return m_timsId;
238}
239void
241 double C0,
242 double C1,
243 double C2,
244 double C3,
245 double C4,
246 [[maybe_unused]] double C5,
247 double C6,
248 double C7,
249 double C8,
250 double C9)
251{
252 if(tims_model_type != 2)
253 {
254 throw pappso::PappsoException(QObject::tr(
255 "ERROR in TimsFrame::setTimsCalibration tims_model_type != 2"));
256 }
257 m_timsDvStart = C2; // C2 from TimsCalibration
258 m_timsTtrans = C4; // C4 from TimsCalibration
259 m_timsNdelay = C0; // C0 from TimsCalibration
260 m_timsVmin = C8; // C8 from TimsCalibration
261 m_timsVmax = C9; // C9 from TimsCalibration
262 m_timsC6 = C6;
263 m_timsC7 = C7;
264
265
267 (C3 - m_timsDvStart) / C1; // //C3 from TimsCalibration // C2 from
268 // TimsCalibration // C1 from TimsCalibration
269}
270double
272{
273 double v = m_timsDvStart +
274 m_timsSlope * ((double)scanNum - m_timsTtrans - m_timsNdelay);
275
276 if(v < m_timsVmin)
277 {
279 QObject::tr("ERROR in TimsFrame::getVoltageTransformation invalid tims "
280 "calibration, v < m_timsVmin"));
281 }
282
283
284 if(v > m_timsVmax)
285 {
287 QObject::tr("ERROR in TimsFrame::getVoltageTransformation invalid tims "
288 "calibration, v > m_timsVmax"));
289 }
290 return v;
291}
292double
293TimsFrameBase::getDriftTime(std::size_t scanNum) const
294{
295 return (m_accumulationTime / (double)m_scanNumber) * ((double)scanNum);
296}
297
298double
300{
301 return 1 / (m_timsC6 + (m_timsC7 / getVoltageTransformation(scanNum)));
302}
303
304
305std::size_t
307{
308 double temp = 1 / one_over_k0;
309 temp = temp - m_timsC6;
310 temp = temp / m_timsC7;
311 temp = 1 / temp;
312 temp = temp - m_timsDvStart;
313 temp = temp / m_timsSlope + m_timsTtrans + m_timsNdelay;
314 return (std::size_t)std::round(temp);
315}
316
317bool
319{
320 if((m_timsDvStart == other.m_timsDvStart) &&
321 (m_timsTtrans == other.m_timsTtrans) &&
322 (m_timsNdelay == other.m_timsNdelay) && (m_timsVmin == other.m_timsVmin) &&
323 (m_timsVmax == other.m_timsVmax) && (m_timsC6 == other.m_timsC6) &&
324 (m_timsC7 == other.m_timsC7) && (m_timsSlope == other.m_timsSlope))
325 {
326 return true;
327 }
328 return false;
329}
330
331
334 std::map<quint32, quint32> &accumulated_scans) const
335{
336 qDebug();
337 // qDebug();
338 // add flanking peaks
339 pappso::Trace local_trace;
340
341 MzCalibrationInterface *mz_calibration_p =
343
344
345 DataPoint element;
346 for(auto &scan_element : accumulated_scans)
347 {
348 // intensity normalization
349 element.y = ((double)scan_element.second) * 100.0 / m_accumulationTime;
350
351 // mz calibration
352 element.x = mz_calibration_p->getMzFromTofIndex(scan_element.first);
353
354 local_trace.push_back(element);
355 }
356 local_trace.sortX();
357
358 qDebug();
359 // qDebug();
360 return local_trace;
361}
362
365 std::map<quint32, quint32> &accumulated_scans) const
366{
367 qDebug();
368 // qDebug();
369 // add flanking peaks
370 std::vector<quint32> keys;
371 transform(begin(accumulated_scans),
372 end(accumulated_scans),
373 back_inserter(keys),
374 [](std::map<quint32, quint32>::value_type const &pair) {
375 return pair.first;
376 });
377 std::sort(keys.begin(), keys.end());
378 pappso::DataPoint data_point_cumul;
379 data_point_cumul.x = 0;
380 data_point_cumul.y = 0;
381
382 pappso::Trace local_trace;
383
384 MzCalibrationInterface *mz_calibration_p =
386
387
388 quint32 last_key = 0;
389
390 for(quint32 key : keys)
391 {
392 if(key == last_key + 1)
393 {
394 // cumulate
395 if(accumulated_scans[key] > accumulated_scans[last_key])
396 {
397 if(data_point_cumul.x == last_key)
398 {
399 // growing peak
400 data_point_cumul.x = key;
401 data_point_cumul.y += accumulated_scans[key];
402 }
403 else
404 {
405 // new peak
406 // flush
407 if(data_point_cumul.y > 0)
408 {
409 // intensity normalization
410 data_point_cumul.y *= 100.0 / m_accumulationTime;
411
412
413 // mz calibration
414 data_point_cumul.x =
415 mz_calibration_p->getMzFromTofIndex(data_point_cumul.x);
416 local_trace.push_back(data_point_cumul);
417 }
418
419 // new point
420 data_point_cumul.x = key;
421 data_point_cumul.y = accumulated_scans[key];
422 }
423 }
424 else
425 {
426 data_point_cumul.y += accumulated_scans[key];
427 }
428 }
429 else
430 {
431 // flush
432 if(data_point_cumul.y > 0)
433 {
434 // intensity normalization
435 data_point_cumul.y *= 100.0 / m_accumulationTime;
436
437
438 // qDebug() << "raw data x=" << data_point_cumul.x;
439 // mz calibration
440 data_point_cumul.x =
441 mz_calibration_p->getMzFromTofIndex(data_point_cumul.x);
442 // qDebug() << "mz=" << data_point_cumul.x;
443 local_trace.push_back(data_point_cumul);
444 }
445
446 // new point
447 data_point_cumul.x = key;
448 data_point_cumul.y = accumulated_scans[key];
449 }
450
451 last_key = key;
452 }
453 // flush
454 if(data_point_cumul.y > 0)
455 {
456 // intensity normalization
457 data_point_cumul.y *= 100.0 / m_accumulationTime;
458
459
460 // mz calibration
461 data_point_cumul.x =
462 mz_calibration_p->getMzFromTofIndex(data_point_cumul.x);
463 local_trace.push_back(data_point_cumul);
464 }
465
466 local_trace.sortX();
467 qDebug();
468 // qDebug();
469 return local_trace;
470}
471
474{
475 if(msp_mzCalibration == nullptr)
476 {
477
479 QObject::tr("ERROR in %1, %2, %3 msp_mzCalibration is null")
480 .arg(__FILE__)
481 .arg(__FUNCTION__)
482 .arg(__LINE__));
483 }
484 return msp_mzCalibration;
485}
486
487void
489 MzCalibrationInterfaceSPtr mzCalibration)
490{
491
492 if(mzCalibration == nullptr)
493 {
494
496 QObject::tr("ERROR in %1, %2, %3 msp_mzCalibration is null")
497 .arg(__FILE__)
498 .arg(__FUNCTION__)
499 .arg(__LINE__));
500 }
501 msp_mzCalibration = mzCalibration;
502}
503
504
505quint32
507{
508 quint32 max_value = 0;
509 for(quint32 i = 0; i < m_scanNumber; i++)
510 {
511 qDebug() << "m_scanNumber=" << m_scanNumber << " i=" << i;
512 std::vector<quint32> index_list = getScanIndexList(i);
513 auto it = std::max_element(index_list.begin(), index_list.end());
514 if(it != index_list.end())
515 {
516 max_value = std::max(max_value, *it);
517 }
518 }
519 return max_value;
520}
521
522std::vector<quint32>
523TimsFrameBase::getScanIndexList(std::size_t scanNum) const
524{
525 throw PappsoException(
526 QObject::tr(
527 "ERROR unable to getScanIndexList in TimsFrameBase for scan number %1")
528 .arg(scanNum));
529}
530
531
532std::vector<quint32>
533TimsFrameBase::getScanIntensities(std::size_t scanNum) const
534{
535 throw PappsoException(
536 QObject::tr(
537 "ERROR unable to getScanIntensities in TimsFrameBase for scan number %1")
538 .arg(scanNum));
539}
540
541Trace
543 std::size_t mz_index_lower_bound,
544 std::size_t mz_index_upper_bound,
545 XicExtractMethod method) const
546{
547 Trace im_trace;
548 DataPoint data_point;
549 for(quint32 i = 0; i < m_scanNumber; i++)
550 {
551 data_point.x = i;
552 data_point.y = 0;
553 qDebug() << "m_scanNumber=" << m_scanNumber << " i=" << i;
554 std::vector<quint32> index_list = getScanIndexList(i);
555 auto it_lower = std::find_if(index_list.begin(),
556 index_list.end(),
557 [mz_index_lower_bound](quint32 to_compare) {
558 if(to_compare < mz_index_lower_bound)
559 {
560 return false;
561 }
562 return true;
563 });
564
565
566 if(it_lower == index_list.end())
567 {
568 }
569 else
570 {
571
572
573 auto it_upper =
574 std::find_if(index_list.begin(),
575 index_list.end(),
576 [mz_index_upper_bound](quint32 to_compare) {
577 if(mz_index_upper_bound >= to_compare)
578 {
579 return false;
580 }
581 return true;
582 });
583 std::vector<quint32> intensity_list = getScanIntensities(i);
584 for(int j = std::distance(index_list.begin(), it_lower);
585 j < std::distance(index_list.begin(), it_upper);
586 j++)
587 {
588 if(method == XicExtractMethod::sum)
589 {
590 data_point.y += intensity_list[j];
591 }
592 else
593 {
594 data_point.y =
595 std::max((double)intensity_list[j], data_point.y);
596 }
597 }
598 }
599 im_trace.push_back(data_point);
600 }
601 qDebug();
602 return im_trace;
603}
604
605std::map<quint32, quint32> &
606TimsFrameBase::downsizeMzRawMap(std::size_t mzindex_merge_window,
607 std::map<quint32, quint32> &rawSpectrum) const
608{
609 std::map<quint32, quint32> new_spectrum;
610
611 for(auto &pair_mz_intensity : rawSpectrum)
612 {
613 quint32 mzkey = (pair_mz_intensity.first / mzindex_merge_window);
614 mzkey = (mzkey * mzindex_merge_window) + (mzindex_merge_window / 2);
615 auto it = new_spectrum.insert({mzkey, pair_mz_intensity.second});
616 if(it.second == false)
617 {
618 it.first->second += pair_mz_intensity.second;
619 }
620 }
621 rawSpectrum = new_spectrum;
622 return rawSpectrum;
623}
624// namespace pappso
625} // 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
virtual quint64 cumulateSingleScanIntensities(std::size_t scanNum) const
double getVoltageTransformation(std::size_t scanNum) const
get voltage for a given scan number
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
virtual Trace cumulateScansToTraceMzDownResolution(std::size_t mzindex_merge_window, std::size_t scanNumBegin, std::size_t scanNumEnd, quint32 &minimum_index, quint32 &maximum_index) const
cumulate spectrum given a scan number range need the binary file The intensities are normalized with ...
std::size_t getScanNumFromOneOverK0(double one_over_k0) const
get the scan number from a given 1/Ko mobility value
virtual Trace cumulateScansToTrace(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 ...
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(SortOrder sort_order=SortOrder::ascending)
Definition trace.cpp:1086
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
XicExtractMethod
Definition types.h:240
@ 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