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