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