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