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