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