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 <QtEndian>
32#include <memory>
33#include <solvers.h>
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 qDebug();
338}
339
340
341Trace
342TimsFrame::cumulateScanToTrace(std::size_t scanNumBegin,
343 std::size_t scanNumEnd) const
344{
345 qDebug();
346
347 Trace new_trace;
348
349 try
350 {
351 if(m_timsDataFrame.size() == 0)
352 return new_trace;
353 std::map<quint32, quint32> raw_spectrum;
354 // double local_accumulationTime = 0;
355
356 std::size_t imax = scanNumEnd + 1;
357 qDebug();
358 for(std::size_t i = scanNumBegin; i < imax; i++)
359 {
360 qDebug() << i;
361 cumulateScan(i, raw_spectrum);
362 qDebug() << i;
363 // local_accumulationTime += m_accumulationTime;
364 }
365 qDebug();
366 pappso::DataPoint data_point_cumul;
367
368
369 MzCalibrationInterface *mz_calibration_p =
371
372
373 for(std::pair<quint32, quint32> pair_tof_intensity : raw_spectrum)
374 {
375 data_point_cumul.x =
376 mz_calibration_p->getMzFromTofIndex(pair_tof_intensity.first);
377 // normalization
378 data_point_cumul.y =
379 pair_tof_intensity.second * ((double)100.0 / m_accumulationTime);
380 new_trace.push_back(data_point_cumul);
381 }
382 new_trace.sortX();
383 qDebug();
384 }
385
386 catch(std::exception &error)
387 {
388 qDebug() << QString(
389 "Failure in TimsFrame::cumulateScanToTrace %1 to %2 :\n %3")
390 .arg(scanNumBegin, scanNumEnd)
391 .arg(error.what());
392 }
393 return new_trace;
394}
395
396
397void
398TimsFrame::cumulateScansInRawMap(std::map<quint32, quint32> &rawSpectrum,
399 std::size_t scanNumBegin,
400 std::size_t scanNumEnd) const
401{
402 qDebug() << "begin scanNumBegin=" << scanNumBegin
403 << " scanNumEnd=" << scanNumEnd;
404
405 if(m_timsDataFrame.size() == 0)
406 return;
407 try
408 {
409
410 std::size_t imax = scanNumEnd + 1;
411 qDebug();
412 for(std::size_t i = scanNumBegin; i < imax; i++)
413 {
414 qDebug() << i;
415 cumulateScan(i, rawSpectrum);
416 qDebug() << i;
417 // local_accumulationTime += m_accumulationTime;
418 }
419 }
420
421 catch(std::exception &error)
422 {
423 qDebug() << QString("Failure in %1 %2 to %3 :\n %4")
424 .arg(__FUNCTION__)
425 .arg(scanNumBegin)
426 .arg(scanNumEnd)
427 .arg(error.what());
428 }
429 qDebug() << "end";
430}
431
432
434TimsFrame::getMassSpectrumCstSPtr(std::size_t scanNum) const
435{
436 qDebug();
437 return getMassSpectrumSPtr(scanNum);
438}
439
441TimsFrame::getMassSpectrumSPtr(std::size_t scanNum) const
442{
443
444 qDebug() << " scanNum=" << scanNum;
445
446 checkScanNum(scanNum);
447
448 qDebug();
449
450 pappso::MassSpectrumSPtr mass_spectrum_sptr =
451 std::make_shared<pappso::MassSpectrum>();
452 // std::vector<DataPoint>
453
454 if(m_timsDataFrame.size() == 0)
455 return mass_spectrum_sptr;
456 qDebug();
457
458 std::size_t size = getNbrPeaks(scanNum);
459
460 std::size_t offset = getScanOffset(scanNum);
461
462 MzCalibrationInterface *mz_calibration_p =
464
465
466 qint32 previous = -1;
467 qint32 tof_index;
468 // std::vector<quint32> index_list;
469 DataPoint data_point;
470 for(std::size_t i = 0; i < size; i++)
471 {
472 tof_index =
473 (*(quint32 *)((m_timsDataFrame.constData() + (offset * 4) + (i * 8))) +
474 previous);
475 data_point.y = (*(quint32 *)(m_timsDataFrame.constData() + (offset * 4) +
476 (i * 8) + 4));
477
478 // intensity normalization
479 data_point.y *= 100.0 / m_accumulationTime;
480
481 previous = tof_index;
482
483
484 // mz calibration
485 data_point.x = mz_calibration_p->getMzFromTofIndex(tof_index);
486 mass_spectrum_sptr.get()->push_back(data_point);
487 }
488 qDebug();
489 return mass_spectrum_sptr;
490}
491
492
493void
495 std::vector<XicCoordTims *>::iterator &itXicListbegin,
496 std::vector<XicCoordTims *>::iterator &itXicListend,
497 XicExtractMethod method) const
498{
499 qDebug() << std::distance(itXicListbegin, itXicListend);
500 std::vector<TimsFrame::XicComputeStructure> tmp_xic_list;
501
502 for(auto it = itXicListbegin; it != itXicListend; it++)
503 {
504 tmp_xic_list.push_back(TimsFrame::XicComputeStructure(this, **it));
505
506 qDebug() << " tmp_xic_struct.mobilityIndexBegin="
507 << tmp_xic_list.back().mobilityIndexBegin
508 << " tmp_xic_struct.mobilityIndexEnd="
509 << tmp_xic_list.back().mobilityIndexEnd;
510
511 qDebug() << " tmp_xic_struct.mzIndexLowerBound="
512 << tmp_xic_list.back().mzIndexLowerBound
513 << " tmp_xic_struct.mzIndexUpperBound="
514 << tmp_xic_list.back().mzIndexUpperBound;
515 }
516 if(tmp_xic_list.size() == 0)
517 return;
518 /*
519 std::sort(tmp_xic_list.begin(), tmp_xic_list.end(), [](const
520 TimsXicStructure &a, const TimsXicStructure &b) { return
521 a.mobilityIndexBegin < b.mobilityIndexBegin;
522 });
523 */
524 std::vector<std::size_t> unique_scan_num_list;
525 for(auto &&struct_xic : tmp_xic_list)
526 {
527 for(std::size_t scan = struct_xic.mobilityIndexBegin;
528 (scan <= struct_xic.mobilityIndexEnd) && (scan < m_scanNumber);
529 scan++)
530 {
531 unique_scan_num_list.push_back(scan);
532 }
533 }
534 std::sort(unique_scan_num_list.begin(), unique_scan_num_list.end());
535 auto it_scan_num_end =
536 std::unique(unique_scan_num_list.begin(), unique_scan_num_list.end());
537 auto it_scan_num = unique_scan_num_list.begin();
538
539 while(it_scan_num != it_scan_num_end)
540 {
541 TraceSPtr ms_spectrum = getRawTraceSPtr(*it_scan_num);
542 // qDebug() << ms_spectrum.get()->toString();
543 for(auto &&tmp_xic_struct : tmp_xic_list)
544 {
545 if(((*it_scan_num) >= tmp_xic_struct.mobilityIndexBegin) &&
546 ((*it_scan_num) <= tmp_xic_struct.mobilityIndexEnd))
547 {
548 if(method == XicExtractMethod::max)
549 {
550 tmp_xic_struct.tmpIntensity +=
551 ms_spectrum.get()->maxY(tmp_xic_struct.mzIndexLowerBound,
552 tmp_xic_struct.mzIndexUpperBound);
553
554 qDebug() << "tmp_xic_struct.tmpIntensity="
555 << tmp_xic_struct.tmpIntensity;
556 }
557 else
558 {
559 // sum
560 tmp_xic_struct.tmpIntensity +=
561 ms_spectrum.get()->sumY(tmp_xic_struct.mzIndexLowerBound,
562 tmp_xic_struct.mzIndexUpperBound);
563 qDebug() << "tmp_xic_struct.tmpIntensity="
564 << tmp_xic_struct.tmpIntensity;
565 }
566 }
567 }
568 it_scan_num++;
569 }
570
571 for(auto &&tmp_xic_struct : tmp_xic_list)
572 {
573 if(tmp_xic_struct.tmpIntensity != 0)
574 {
575 qDebug() << tmp_xic_struct.xic_ptr;
576 tmp_xic_struct.xic_ptr->push_back(
577 {m_time, tmp_xic_struct.tmpIntensity});
578 }
579 }
580
581 qDebug();
582}
583
584
586TimsFrame::getRawTraceSPtr(std::size_t scanNum) const
587{
588
589 // qDebug();
590
591 pappso::TraceSPtr trace_sptr = std::make_shared<pappso::Trace>();
592 // std::vector<DataPoint>
593
594 if(m_timsDataFrame.size() == 0)
595 return trace_sptr;
596 // qDebug();
597
598 std::size_t size = getNbrPeaks(scanNum);
599
600 std::size_t offset = getScanOffset(scanNum);
601
602 qint32 previous = -1;
603 std::vector<quint32> index_list;
604 for(std::size_t i = 0; i < size; i++)
605 {
606 DataPoint data_point(
607 (*(quint32 *)((m_timsDataFrame.constData() + (offset * 4) + (i * 8))) +
608 previous),
609 (*(quint32 *)(m_timsDataFrame.constData() + (offset * 4) + (i * 8) +
610 4)));
611
612 // intensity normalization
613 data_point.y *= 100.0 / m_accumulationTime;
614
615 previous = data_point.x;
616 trace_sptr.get()->push_back(data_point);
617 }
618 // qDebug();
619 return trace_sptr;
620}
621
622} // 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:441
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:398
virtual Trace cumulateScanToTrace(std::size_t scanNumBegin, std::size_t scanNumEnd) const override
cumulate scan list into a trace
Definition: timsframe.cpp:342
void extractTimsXicListInRtRange(std::vector< XicCoordTims * >::iterator &itXicListbegin, std::vector< XicCoordTims * >::iterator &itXicListend, XicExtractMethod method) const
Definition: timsframe.cpp:494
virtual pappso::MassSpectrumCstSPtr getMassSpectrumCstSPtr(std::size_t scanNum) const
get the mass spectrum corresponding to a scan number
Definition: timsframe.cpp:434
virtual pappso::TraceSPtr getRawTraceSPtr(std::size_t scanNum) const
get the raw index tof_index and intensities (normalized)
Definition: timsframe.cpp:586
A simple container of DataPoint instances.
Definition: trace.h:38
Q_INVOKABLE void sortX()
Definition: trace.cpp:983
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:25
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: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: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