libpappsomspp
Library for mass spectrometry
timsframe.cpp
Go to the documentation of this file.
1/**
2 * \file pappsomspp/vendors/tims/timsframe.cpp
3 * \date 23/08/2019
4 * \author Olivier Langella
5 * \brief handle a single Bruker's TimsTof frame
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
28#include "timsframe.h"
29#include "../../../pappsomspp/pappsoexception.h"
30#include <QDebug>
31#include <QObject>
32#include <QtEndian>
33#include <memory>
34
35
36namespace pappso
37{
38
40 const TimsFrame *fram_p, const XicCoordTims &xic_struct)
41{
42 xic_ptr = xic_struct.xicSptr.get();
43
45 mobilityIndexEnd = xic_struct.scanNumEnd;
47 fram_p->getMzCalibrationInterfaceSPtr().get()->getTofIndexFromMz(
48 xic_struct.mzRange.lower()); // convert mz to raw digitizer value
50 fram_p->getMzCalibrationInterfaceSPtr().get()->getTofIndexFromMz(
51 xic_struct.mzRange.upper());
52 tmpIntensity = 0;
53}
54
55TimsFrame::TimsFrame(std::size_t timsId, quint32 scanNum)
56 : TimsFrameBase(timsId, scanNum)
57{
58 // m_timsDataFrame.resize(10);
59}
60
61TimsFrame::TimsFrame(std::size_t timsId,
62 quint32 scanNum,
63 char *p_bytes,
64 std::size_t len)
65 : TimsFrameBase(timsId, scanNum)
66{
67 // langella@themis:~/developpement/git/bruker/cbuild$
68 // ./src/sample/timsdataSamplePappso
69 // /gorgone/pappso/fichiers_fabricants/Bruker/Demo_TimsTOF_juin2019/Samples/1922001/1922001-1_S-415_Pep_Pur-1ul_Slot1-10_1_2088.d/
70 qDebug() << timsId;
71
72 m_timsDataFrame.resize(len);
73
74 if(p_bytes != nullptr)
75 {
76 unshufflePacket(p_bytes);
77 }
78 else
79 {
80 if(m_scanNumber == 0)
81 {
82
84 QObject::tr("TimsFrame::TimsFrame(%1,%2,nullptr,%3) FAILED")
85 .arg(m_timsId)
86 .arg(m_scanNumber)
87 .arg(len));
88 }
89 }
90}
91
93{
94}
95
97{
98}
99
100
101void
103{
104 qDebug();
105 quint64 len = m_timsDataFrame.size();
106 if(len % 4 != 0)
107 {
109 QObject::tr("TimsFrame::unshufflePacket error: len % 4 != 0"));
110 }
111
112 quint64 nb_uint4 = len / 4;
113
114 char *dest = m_timsDataFrame.data();
115 quint64 src_offset = 0;
116
117 for(quint64 j = 0; j < 4; j++)
118 {
119 for(quint64 i = 0; i < nb_uint4; i++)
120 {
121 dest[(i * 4) + j] = src[src_offset];
122 src_offset++;
123 }
124 }
125 qDebug();
126}
127
128std::size_t
129TimsFrame::getNbrPeaks(std::size_t scanNum) const
130{
131 if(m_timsDataFrame.size() == 0)
132 return 0;
133 /*
134 if(scanNum == 0)
135 {
136 quint32 res = (*(quint32 *)(m_timsDataFrame.constData() + 4)) -
137 (*(quint32 *)(m_timsDataFrame.constData()-4));
138 return res / 2;
139 }*/
140 if(scanNum == (m_scanNumber - 1))
141 {
142 auto nb_uint4 = m_timsDataFrame.size() / 4;
143
144 std::size_t cumul = 0;
145 for(quint32 i = 0; i < m_scanNumber; i++)
146 {
147 cumul += (*(quint32 *)(m_timsDataFrame.constData() + (i * 4)));
148 }
149 return (nb_uint4 - cumul) / 2;
150 }
151 checkScanNum(scanNum);
152
153 // quint32 *res = (quint32 *)(m_timsDataFrame.constData() + (scanNum * 4));
154 // qDebug() << " res=" << *res;
155 return (*(quint32 *)(m_timsDataFrame.constData() + ((scanNum + 1) * 4))) / 2;
156}
157
158std::size_t
159TimsFrame::getScanOffset(std::size_t scanNum) const
160{
161 std::size_t offset = 0;
162 for(std::size_t i = 0; i < (scanNum + 1); i++)
163 {
164 offset += (*(quint32 *)(m_timsDataFrame.constData() + (i * 4)));
165 }
166 return offset;
167}
168
169
170std::vector<quint32>
171TimsFrame::getScanIndexList(std::size_t scanNum) const
172{
173 qDebug();
174 checkScanNum(scanNum);
175 std::vector<quint32> scan_tof;
176
177 if(m_timsDataFrame.size() == 0)
178 return scan_tof;
179 scan_tof.resize(getNbrPeaks(scanNum));
180
181 std::size_t offset = getScanOffset(scanNum);
182
183 qint32 previous = -1;
184 for(std::size_t i = 0; i < scan_tof.size(); i++)
185 {
186 scan_tof[i] =
187 (*(quint32 *)(m_timsDataFrame.constData() + (offset * 4) + (i * 8))) +
188 previous;
189 previous = scan_tof[i];
190 }
191 qDebug();
192 return scan_tof;
193}
194
195std::vector<quint32>
196TimsFrame::getScanIntensities(std::size_t scanNum) const
197{
198 qDebug();
199 checkScanNum(scanNum);
200 std::vector<quint32> scan_intensities;
201
202 if(m_timsDataFrame.size() == 0)
203 return scan_intensities;
204
205 scan_intensities.resize(getNbrPeaks(scanNum));
206
207 std::size_t offset = getScanOffset(scanNum);
208
209 for(std::size_t i = 0; i < scan_intensities.size(); i++)
210 {
211 scan_intensities[i] = (*(quint32 *)(m_timsDataFrame.constData() +
212 (offset * 4) + (i * 8) + 4));
213 }
214 qDebug();
215 return scan_intensities;
216}
217
218
219quint64
221{
222 qDebug();
223
224 quint64 summed_intensities = 0;
225
226 if(m_timsDataFrame.size() == 0)
227 return summed_intensities;
228 // checkScanNum(scanNum);
229
230 std::size_t size = getNbrPeaks(scanNum);
231
232 std::size_t offset = getScanOffset(scanNum);
233
234 qint32 previous = -1;
235
236 for(std::size_t i = 0; i < size; i++)
237 {
238 quint32 x =
239 (*(quint32 *)((m_timsDataFrame.constData() + (offset * 4) + (i * 8))) +
240 previous);
241
242 quint32 y = (*(quint32 *)(m_timsDataFrame.constData() + (offset * 4) +
243 (i * 8) + 4));
244
245 previous = x;
246
247 summed_intensities += y;
248 }
249
250 // Normalization over the accumulation time for this frame.
251 summed_intensities *= ((double)100.0 / m_accumulationTime);
252
253 qDebug();
254
255 return summed_intensities;
256}
257
258
259/**
260 * @brief ...
261 *
262 * @param scanNumBegin p_scanNumBegin:...
263 * @param scanNumEnd p_scanNumEnd:...
264 * @return quint64
265 */
266quint64
267TimsFrame::cumulateScansIntensities(std::size_t scanNumBegin,
268 std::size_t scanNumEnd) const
269{
270 quint64 summed_intensities = 0;
271
272 // qDebug() << "begin scanNumBegin =" << scanNumBegin
273 //<< "scanNumEnd =" << scanNumEnd;
274
275 if(m_timsDataFrame.size() == 0)
276 return summed_intensities;
277
278 try
279 {
280 std::size_t imax = scanNumEnd + 1;
281
282 for(std::size_t i = scanNumBegin; i < imax; i++)
283 {
284 qDebug() << i;
285 summed_intensities += cumulateSingleScanIntensities(i);
286 qDebug() << i;
287 }
288 }
289 catch(std::exception &error)
290 {
291 qDebug() << QString("Failure in %1 %2 to %3 :\n %4")
292 .arg(__FUNCTION__)
293 .arg(scanNumBegin)
294 .arg(scanNumEnd)
295 .arg(error.what());
296 }
297
298 // qDebug() << "end";
299
300 return summed_intensities;
301}
302
303
304void
305TimsFrame::cumulateScan(std::size_t scanNum,
306 std::map<quint32, quint32> &accumulate_into) const
307{
308 // qDebug();
309
310 if(m_timsDataFrame.size() == 0)
311 return;
312 // checkScanNum(scanNum);
313
314 std::size_t size = getNbrPeaks(scanNum);
315
316 std::size_t offset = getScanOffset(scanNum);
317
318 qint32 previous = -1;
319 for(std::size_t i = 0; i < size; i++)
320 {
321 quint32 x =
322 (*(quint32 *)((m_timsDataFrame.constData() + (offset * 4) + (i * 8))) +
323 previous);
324 quint32 y = (*(quint32 *)(m_timsDataFrame.constData() + (offset * 4) +
325 (i * 8) + 4));
326
327 previous = x;
328
329 auto ret = accumulate_into.insert(std::pair<quint32, quint32>(x, y));
330
331 if(ret.second == false)
332 {
333 // already existed : cumulate
334 ret.first->second += y;
335 }
336 }
337
338 // qDebug();
339}
340
341
342Trace
343TimsFrame::cumulateScanToTrace(std::size_t scanNumBegin,
344 std::size_t scanNumEnd) const
345{
346 // qDebug();
347
348 Trace new_trace;
349
350 try
351 {
352 if(m_timsDataFrame.size() == 0)
353 return new_trace;
354 std::map<quint32, quint32> raw_spectrum;
355 // double local_accumulationTime = 0;
356
357 std::size_t imax = scanNumEnd + 1;
358 qDebug();
359 for(std::size_t i = scanNumBegin; i < imax; i++)
360 {
361 // qDebug() << i;
362 cumulateScan(i, raw_spectrum);
363 // qDebug() << i;
364
365 // local_accumulationTime += m_accumulationTime;
366 }
367
368 // qDebug();
369
370 pappso::DataPoint data_point_cumul;
371
372
373 MzCalibrationInterface *mz_calibration_p =
375
376
377 for(std::pair<quint32, quint32> pair_tof_intensity : raw_spectrum)
378 {
379 data_point_cumul.x =
380 mz_calibration_p->getMzFromTofIndex(pair_tof_intensity.first);
381 // normalization
382 data_point_cumul.y =
383 pair_tof_intensity.second * ((double)100.0 / m_accumulationTime);
384 new_trace.push_back(data_point_cumul);
385 }
386 new_trace.sortX();
387
388 // qDebug();
389 }
390
391 catch(std::exception &error)
392 {
393 qDebug() << QString(
394 "Failure in TimsFrame::cumulateScanToTrace %1 to %2 :\n %3")
395 .arg(scanNumBegin, scanNumEnd)
396 .arg(error.what());
397 }
398 return new_trace;
399}
400
401
402void
403TimsFrame::cumulateScansInRawMap(std::map<quint32, quint32> &rawSpectrum,
404 std::size_t scanNumBegin,
405 std::size_t scanNumEnd) const
406{
407 // qDebug() << "begin scanNumBegin=" << scanNumBegin
408 //<< " scanNumEnd=" << scanNumEnd;
409
410 if(m_timsDataFrame.size() == 0)
411 return;
412 try
413 {
414
415 std::size_t imax = scanNumEnd + 1;
416 qDebug();
417 for(std::size_t i = scanNumBegin; i < imax; i++)
418 {
419 qDebug() << i;
420 cumulateScan(i, rawSpectrum);
421 qDebug() << i;
422 // local_accumulationTime += m_accumulationTime;
423 }
424 }
425
426 catch(std::exception &error)
427 {
428 qDebug() << QString("Failure in %1 %2 to %3 :\n %4")
429 .arg(__FUNCTION__)
430 .arg(scanNumBegin)
431 .arg(scanNumEnd)
432 .arg(error.what());
433 }
434
435 // qDebug() << "end";
436}
437
438
440TimsFrame::getMassSpectrumCstSPtr(std::size_t scanNum) const
441{
442 // qDebug();
443
444 return getMassSpectrumSPtr(scanNum);
445}
446
448TimsFrame::getMassSpectrumSPtr(std::size_t scanNum) const
449{
450
451 // qDebug() << " scanNum=" << scanNum;
452
453 checkScanNum(scanNum);
454
455 // qDebug();
456
457 pappso::MassSpectrumSPtr mass_spectrum_sptr =
458 std::make_shared<pappso::MassSpectrum>();
459 // std::vector<DataPoint>
460
461 if(m_timsDataFrame.size() == 0)
462 return mass_spectrum_sptr;
463
464 // qDebug();
465
466 std::size_t size = getNbrPeaks(scanNum);
467
468 std::size_t offset = getScanOffset(scanNum);
469
470 MzCalibrationInterface *mz_calibration_p =
472
473
474 qint32 previous = -1;
475 qint32 tof_index;
476 // std::vector<quint32> index_list;
477 DataPoint data_point;
478 for(std::size_t i = 0; i < size; i++)
479 {
480 tof_index =
481 (*(quint32 *)((m_timsDataFrame.constData() + (offset * 4) + (i * 8))) +
482 previous);
483 data_point.y = (*(quint32 *)(m_timsDataFrame.constData() + (offset * 4) +
484 (i * 8) + 4));
485
486 // intensity normalization
487 data_point.y *= 100.0 / m_accumulationTime;
488
489 previous = tof_index;
490
491
492 // mz calibration
493 data_point.x = mz_calibration_p->getMzFromTofIndex(tof_index);
494 mass_spectrum_sptr.get()->push_back(data_point);
495 }
496
497 // qDebug();
498
499 return mass_spectrum_sptr;
500}
501
502
503void
505 std::vector<XicCoordTims *>::iterator &itXicListbegin,
506 std::vector<XicCoordTims *>::iterator &itXicListend,
507 XicExtractMethod method) const
508{
509 qDebug() << std::distance(itXicListbegin, itXicListend);
510
511 std::vector<TimsFrame::XicComputeStructure> tmp_xic_list;
512
513 for(auto it = itXicListbegin; it != itXicListend; it++)
514 {
515 tmp_xic_list.push_back(TimsFrame::XicComputeStructure(this, **it));
516
517 qDebug() << " tmp_xic_struct.mobilityIndexBegin="
518 << tmp_xic_list.back().mobilityIndexBegin
519 << " tmp_xic_struct.mobilityIndexEnd="
520 << tmp_xic_list.back().mobilityIndexEnd;
521
522 qDebug() << " tmp_xic_struct.mzIndexLowerBound="
523 << tmp_xic_list.back().mzIndexLowerBound
524 << " tmp_xic_struct.mzIndexUpperBound="
525 << tmp_xic_list.back().mzIndexUpperBound;
526 }
527 if(tmp_xic_list.size() == 0)
528 return;
529 /*
530 std::sort(tmp_xic_list.begin(), tmp_xic_list.end(), [](const
531 TimsXicStructure &a, const TimsXicStructure &b) { return
532 a.mobilityIndexBegin < b.mobilityIndexBegin;
533 });
534 */
535 std::vector<std::size_t> unique_scan_num_list;
536 for(auto &&struct_xic : tmp_xic_list)
537 {
538 for(std::size_t scan = struct_xic.mobilityIndexBegin;
539 (scan <= struct_xic.mobilityIndexEnd) && (scan < m_scanNumber);
540 scan++)
541 {
542 unique_scan_num_list.push_back(scan);
543 }
544 }
545 std::sort(unique_scan_num_list.begin(), unique_scan_num_list.end());
546 auto it_scan_num_end =
547 std::unique(unique_scan_num_list.begin(), unique_scan_num_list.end());
548 auto it_scan_num = unique_scan_num_list.begin();
549
550 while(it_scan_num != it_scan_num_end)
551 {
552 TraceSPtr ms_spectrum = getRawTraceSPtr(*it_scan_num);
553 // qDebug() << ms_spectrum.get()->toString();
554 for(auto &&tmp_xic_struct : tmp_xic_list)
555 {
556 if(((*it_scan_num) >= tmp_xic_struct.mobilityIndexBegin) &&
557 ((*it_scan_num) <= tmp_xic_struct.mobilityIndexEnd))
558 {
559 if(method == XicExtractMethod::max)
560 {
561 tmp_xic_struct.tmpIntensity +=
562 ms_spectrum.get()->maxY(tmp_xic_struct.mzIndexLowerBound,
563 tmp_xic_struct.mzIndexUpperBound);
564
565 qDebug() << "tmp_xic_struct.tmpIntensity="
566 << tmp_xic_struct.tmpIntensity;
567 }
568 else
569 {
570 // sum
571 tmp_xic_struct.tmpIntensity +=
572 ms_spectrum.get()->sumY(tmp_xic_struct.mzIndexLowerBound,
573 tmp_xic_struct.mzIndexUpperBound);
574 qDebug() << "tmp_xic_struct.tmpIntensity="
575 << tmp_xic_struct.tmpIntensity;
576 }
577 }
578 }
579 it_scan_num++;
580 }
581
582 for(auto &&tmp_xic_struct : tmp_xic_list)
583 {
584 if(tmp_xic_struct.tmpIntensity != 0)
585 {
586 qDebug() << tmp_xic_struct.xic_ptr;
587 tmp_xic_struct.xic_ptr->push_back(
588 {m_time, tmp_xic_struct.tmpIntensity});
589 }
590 }
591
592 qDebug();
593}
594
595
597TimsFrame::getRawTraceSPtr(std::size_t scanNum) const
598{
599
600 // qDebug();
601
602 pappso::TraceSPtr trace_sptr = std::make_shared<pappso::Trace>();
603 // std::vector<DataPoint>
604
605 if(m_timsDataFrame.size() == 0)
606 return trace_sptr;
607 // qDebug();
608
609 std::size_t size = getNbrPeaks(scanNum);
610
611 std::size_t offset = getScanOffset(scanNum);
612
613 qint32 previous = -1;
614 std::vector<quint32> index_list;
615 for(std::size_t i = 0; i < size; i++)
616 {
617 DataPoint data_point(
618 (*(quint32 *)((m_timsDataFrame.constData() + (offset * 4) + (i * 8))) +
619 previous),
620 (*(quint32 *)(m_timsDataFrame.constData() + (offset * 4) + (i * 8) +
621 4)));
622
623 // intensity normalization
624 data_point.y *= 100.0 / m_accumulationTime;
625
626 previous = data_point.x;
627 trace_sptr.get()->push_back(data_point);
628 }
629 // qDebug();
630 return trace_sptr;
631}
632
633} // namespace pappso
virtual double getMzFromTofIndex(quint32 tof_index)=0
get m/z from time of flight raw index
pappso_double lower() const
Definition: mzrange.h:71
pappso_double upper() const
Definition: mzrange.h:77
double m_accumulationTime
accumulation time in milliseconds
double m_time
retention time
quint32 m_scanNumber
total number of scans contained in this frame
std::size_t m_timsId
Tims frame database id (the SQL identifier of this frame)
virtual const MzCalibrationInterfaceSPtr & getMzCalibrationInterfaceSPtr() const final
get the MzCalibration model to compute mz and TOF for this frame
bool checkScanNum(std::size_t scanNum) const
check that this scan number exists
QByteArray m_timsDataFrame
Definition: timsframe.h:187
virtual std::vector< quint32 > getScanIndexList(std::size_t scanNum) const override
get raw index list for one given scan index are not TOF nor m/z, just index on digitizer
Definition: timsframe.cpp:171
virtual std::size_t getNbrPeaks(std::size_t scanNum) const override
get the number of peaks in this spectrum need the binary file
Definition: timsframe.cpp:129
virtual ~TimsFrame()
Definition: timsframe.cpp:96
virtual quint64 cumulateScansIntensities(std::size_t scanNumBegin, std::size_t scanNumEnd) const override
...
Definition: timsframe.cpp:267
virtual void cumulateScan(std::size_t scanNum, std::map< quint32, quint32 > &accumulate_into) const
cumulate a scan into a map
Definition: timsframe.cpp:305
TimsFrame(std::size_t timsId, quint32 scanNum, char *p_bytes, std::size_t len)
Definition: timsframe.cpp:61
virtual pappso::MassSpectrumSPtr getMassSpectrumSPtr(std::size_t scanNum) const override
get Mass spectrum with peaks for this scan number need the binary file
Definition: timsframe.cpp:448
virtual quint64 cumulateSingleScanIntensities(std::size_t scanNum) const override
Definition: timsframe.cpp:220
virtual std::vector< quint32 > getScanIntensities(std::size_t scanNum) const override
get raw intensities without transformation from one scan it needs intensity normalization
Definition: timsframe.cpp:196
void unshufflePacket(const char *src)
unshuffle data packet of tims compression type 2
Definition: timsframe.cpp:102
std::size_t getScanOffset(std::size_t scanNum) const
get offset for this spectrum in the binary file
Definition: timsframe.cpp:159
virtual void cumulateScansInRawMap(std::map< quint32, quint32 > &rawSpectrum, std::size_t scanNumBegin, std::size_t scanNumEnd) const override
cumulate scan list into a trace into a raw spectrum map
Definition: timsframe.cpp:403
virtual Trace cumulateScanToTrace(std::size_t scanNumBegin, std::size_t scanNumEnd) const override
cumulate scan list into a trace
Definition: timsframe.cpp:343
void extractTimsXicListInRtRange(std::vector< XicCoordTims * >::iterator &itXicListbegin, std::vector< XicCoordTims * >::iterator &itXicListend, XicExtractMethod method) const
Definition: timsframe.cpp:504
virtual pappso::MassSpectrumCstSPtr getMassSpectrumCstSPtr(std::size_t scanNum) const
get the mass spectrum corresponding to a scan number
Definition: timsframe.cpp:440
virtual pappso::TraceSPtr getRawTraceSPtr(std::size_t scanNum) const
get the raw index tof_index and intensities (normalized)
Definition: timsframe.cpp:597
A simple container of DataPoint instances.
Definition: trace.h:148
void sortX(SortOrder sort_order=SortOrder::ascending)
Definition: trace.cpp:1086
tries to keep as much as possible monoisotopes, removing any possible C13 peaks and changes multichar...
Definition: aa.cpp:39
std::shared_ptr< Trace > TraceSPtr
Definition: trace.h:135
std::shared_ptr< const MassSpectrum > MassSpectrumCstSPtr
Definition: massspectrum.h:55
std::shared_ptr< MassSpectrum > MassSpectrumSPtr
Definition: massspectrum.h:54
XicExtractMethod
Definition: types.h:240
@ max
maximum of intensities
pappso_double x
Definition: datapoint.h:23
pappso_double y
Definition: datapoint.h:24
XicComputeStructure(const TimsFrame *fram_p, const XicCoordTims &xic_struct)
Definition: timsframe.cpp:39
coordinates of the XIC to extract and the resulting XIC after extraction
Definition: xiccoordtims.h:51
std::size_t scanNumEnd
mobility index end
Definition: xiccoordtims.h:98
std::size_t scanNumBegin
mobility index begin
Definition: xiccoordtims.h:94
XicSPtr xicSptr
extracted xic
Definition: xiccoord.h:130
MzRange mzRange
the mass to extract
Definition: xiccoord.h:120
handle a single Bruker's TimsTof frame