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 XicCoordTims &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, quint32 scanNum)
58  : TimsFrameBase(timsId, scanNum)
59 {
60  // m_timsDataFrame.resize(10);
61 }
62 
63 TimsFrame::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 
103 void
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 
130 std::size_t
131 TimsFrame::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 
160 std::size_t
161 TimsFrame::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 
172 std::vector<quint32>
173 TimsFrame::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 
197 std::vector<quint32>
198 TimsFrame::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 
221 void
222 TimsFrame::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 
259 Trace
260 TimsFrame::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 
315 void
316 TimsFrame::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 
352 TimsFrame::getMassSpectrumCstSPtr(std::size_t scanNum) const
353 {
354  qDebug();
355  return getMassSpectrumSPtr(scanNum);
356 }
357 
359 TimsFrame::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 
411 void
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 
504 TimsFrame::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