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