libpappsomspp
Library for mass spectrometry
Loading...
Searching...
No Matches
timsframesmsrunreader.cpp
Go to the documentation of this file.
1/**
2 * \file pappsomspp/msrun/private/timsmsrunreader.h
3 * \date 05/09/2019
4 * \author Olivier Langella
5 * \brief MSrun file reader for native Bruker TimsTOF raw data
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
33#include <QDebug>
34
35using namespace pappso;
36
41
46
47
50 [[maybe_unused]] std::size_t spectrum_index)
51{
53 QObject::tr("Not yet implemented in TimsFramesMsRunReader %1.\n")
54 .arg(__LINE__));
55
57}
58
59
62{
64}
65
66
69 bool want_binary_data) const
70{
71
72 QualifiedMassSpectrum mass_spectrum;
73
75 getMsRunId(), mass_spectrum, spectrum_index, want_binary_data);
76 return mass_spectrum;
77}
78
79
80void
83{
84 // qDebug() << "Reading the spectrum collection with no specific
85 // configuration.";
86 MsRunReadConfig config;
87 readSpectrumCollection2(config, handler);
88}
89
90
91void
94{
95 // qDebug().noquote() << "Reading the spectrum collection with this "
96 // "specific configuration:"
97 // << config.toString();
98
99 std::vector<std::size_t> subset_of_tims_frame_ids;
100
101
102 bool asked_ion_mobility_scan_num_range = false;
103
104 quint32 mobility_scan_num_range_begin = std::numeric_limits<quint32>::max();
105 quint32 mobility_scan_num_range_end = std::numeric_limits<quint32>::max();
106 quint32 mobility_scan_num_range_width = std::numeric_limits<quint32>::max();
107
108 double mobility_one_over_k0 = std::numeric_limits<double>::max();
109 double mobility_one_over_k0_range_begin = std::numeric_limits<double>::max();
110 double mobility_one_over_k0_range_end = std::numeric_limits<double>::max();
111
112 if(!config
113 .getParameterValue(
115 .isNull() &&
116 !config
117 .getParameterValue(
119 .isNull())
120 {
121 mobility_scan_num_range_begin =
122 config
125 .toUInt();
126 mobility_scan_num_range_end =
127 config
130 .toUInt();
131
132 // We need the range width below.
133 mobility_scan_num_range_width =
134 mobility_scan_num_range_end + 1 - mobility_scan_num_range_begin;
135
136 asked_ion_mobility_scan_num_range = true;
137
138 // Be sure to check in the frames loop below that the user might
139 // have asked for an ion mobility range but on the basis of the 1/K0 unit.
140 }
141
142 const std::vector<FrameIdDescr> &frame_id_descr_list =
144
145 // Just for the feedback to the user.
146 std::size_t scan_count = 0;
147
148 for(auto const &frame_record : msp_timsData->getTimsFrameRecordList())
149 {
150 if(handler.shouldStop())
151 {
152 // qDebug() << "The operation was cancelled. Breaking the loop.";
154 QObject::tr("Reading timsTOF data cancelled by the user."));
155 }
156
157 if(frame_record.frame_id == 0)
158 continue;
159
160 if(!config.acceptRetentionTimeInSeconds(frame_record.frame_time))
161 continue;
162
163 std::size_t ms_level = 2;
164 if(frame_record.msms_type == 0)
165 ms_level = 1;
166
167 if(!config.acceptMsLevel(ms_level))
168 continue;
169
170 // FIXME: this might be the place where to actually get the list of
171 // precursors' m/z and charge values.
172
173 subset_of_tims_frame_ids.push_back(frame_record.frame_id);
174
175 if(mobility_scan_num_range_width != std::numeric_limits<int>::max())
176 {
177 scan_count += mobility_scan_num_range_width;
178 }
179 else
180 {
181 scan_count += frame_id_descr_list[frame_record.frame_id].m_scanCount;
182 }
183 }
184
185 // At this point, we have a subset of frame records.
186 std::size_t frame_count = subset_of_tims_frame_ids.size();
187 qDebug() << "The number of retained RT range- and MS level-matching frames : "
188 << frame_count;
189
190 // Inform the handler of the spectrum list so that it can handle feedback to
191 // the user.
192
193 // FIXME:
194 // Either we document the number of frames (because we assume we will
195 // flatten them all)...
196 handler.spectrumListHasSize(frame_count);
197 // Or we document the number of actual scans because we might not flatten
198 // all the frames.
199 // handler.spectrumListHasSize(scan_count);
200
201 // Check for m/z range selection
202 bool asked_mz_range = false;
203 double mz_range_begin = -1;
204 double mz_range_end = -1;
205
207 {
208 asked_mz_range = true;
209
210 mz_range_begin =
212 .toDouble();
213 mz_range_end =
215 .toDouble();
216
217 // qDebug() << "The m/z range asked is: " << mz_range_begin << "--"
218 // << mz_range_end;
219 }
220
221 // Check for m/z resolution downgrading
222 // The idea is that we merge a number of mz indices into a single index,
223 // which is essentially an increase of the m/z bin size, and therefore
224 // of the resolution/definition of the mass spectrum.
225 std::size_t mz_index_merge_window = 0;
226 if(!config
227 .getParameterValue(
229 .isNull())
230 {
231 mz_index_merge_window =
232 config
235 .toUInt();
236
237 // qDebug() << "mz_index_merge_window=" << mz_index_merge_window;
238 }
239
240 // The scan index is the index of the scan in the *whole* mass data file, it
241 // is a sequential number of scans over all the frames.
242 std::size_t scan_index = 0; // iterate in each spectrum
243
244 for(std::size_t tims_frame_id : subset_of_tims_frame_ids)
245 {
246 if(handler.shouldStop())
247 {
248 // qDebug() << "The operation was cancelled. Breaking the loop.";
250 QObject::tr("Reading timsTOF data cancelled by the user."));
251 }
252
253 // qDebug() << "tims_frame_id=" << tims_frame_id;
254
255 const FrameIdDescr &current_frame_record =
256 frame_id_descr_list[tims_frame_id];
257 // qDebug() << "tims_frame_id=" << tims_frame_id;
258
259 // FIXME: from an outsider point of view, cumulated size does not
260 // convey the notion of sequential scan number.
261
262 scan_index = current_frame_record.m_globalScanIndex;
263
264 TimsFrameCstSPtr tims_frame_csp =
266
267 // If the user wants to select 1/Ko values in a given range, we need to
268 // compute the ion mobility scan value starting from that 1/Ko value in
269 // *each* frame. Note that the computed mobility_scan_num_begin and
270 // mobility_scan_num_end would override thoses possibly set with
271 // TimsFramesMsRunReader_mobility_index_begin/end above.
272
273 if(!config
274 .getParameterValue(
276 .isNull() &&
277 !config
278 .getParameterValue(
280 .isNull())
281 {
282 mobility_one_over_k0_range_begin =
283 config
286 .toDouble();
287
288 mobility_one_over_k0_range_end =
289 config
292 .toDouble();
293
294 mobility_scan_num_range_begin =
295 tims_frame_csp.get()->getScanIndexFromOneOverK0(
296 mobility_one_over_k0_range_begin);
297
298 mobility_scan_num_range_end =
299 tims_frame_csp.get()->getScanIndexFromOneOverK0(
300 mobility_one_over_k0_range_end);
301
302 asked_ion_mobility_scan_num_range = true;
303 }
304
305 // Now that we know if the user has asked for an ion mobility range,
306 // either using scan indices or 1/K0 values, we need to double check the
307 // range borders.
308 quint32 count_of_mobility_scans = tims_frame_csp->getTotalNumberOfScans();
309
310 if(asked_ion_mobility_scan_num_range)
311 {
312 if(mobility_scan_num_range_end > (count_of_mobility_scans - 1))
313 {
314 mobility_scan_num_range_end = count_of_mobility_scans - 1;
315 }
316 }
317 else
318 {
319 mobility_scan_num_range_begin = 0;
320 mobility_scan_num_range_end = count_of_mobility_scans - 1;
321 }
322
323 // Now that we know the mobility index range, if we did not set the
324 // mobility one over K0 because that was not the unit used by
325 // the caller, then we can compute these values and set them
326 // later to the qualified mass spectrum parameters.
327 if(mobility_one_over_k0_range_begin == std::numeric_limits<double>::max())
328 mobility_one_over_k0_range_begin =
329 tims_frame_csp->getOneOverK0Transformation(
330 mobility_scan_num_range_begin);
331 if(mobility_one_over_k0_range_end == std::numeric_limits<double>::max())
332 mobility_one_over_k0_range_end =
333 tims_frame_csp->getOneOverK0Transformation(
334 mobility_scan_num_range_end);
335
336 mobility_scan_num_range_width =
337 mobility_scan_num_range_end + 1 - mobility_scan_num_range_begin;
338
339 // We want to provide the inverse mobility for the scan that sits in the
340 // middle of the defined range or the whole range if none is defined..
341 mobility_one_over_k0 = tims_frame_csp.get()->getScanIndexFromOneOverK0(
342 mobility_scan_num_range_begin + (mobility_scan_num_range_width / 2));
343
344 // Now, with or without the peak list, we have to craft a qualified mass
345 // spectrum that will hold all the data about the data in it.
346 QualifiedMassSpectrum qualified_mass_spectrum;
347
348 MassSpectrumId spectrum_id;
349
350 spectrum_id.setSpectrumIndex(tims_frame_id);
351 spectrum_id.setMsRunId(getMsRunId());
352
353 // Can be modified to add bits that might help our case
354 spectrum_id.setNativeId(
355 QString("frame_id=%1 global_scan_index=%2 im_scan_range_begin=%3 "
356 "im_scan_range_end=%4")
357 .arg(tims_frame_id)
358 .arg(scan_index)
359 .arg(mobility_scan_num_range_begin)
360 .arg(mobility_scan_num_range_end));
361
362 qualified_mass_spectrum.setMassSpectrumId(spectrum_id);
363
364 // We want to document the retention time!
365
366 qualified_mass_spectrum.setRtInSeconds(
367 tims_frame_csp.get()->getRtInSeconds());
368
369 // We do want to document the ms level of the spectrum and possibly
370 // the precursor's m/z and charge.
371 unsigned int frame_ms_level = tims_frame_csp.get()->getMsLevel();
372 qualified_mass_spectrum.setMsLevel(frame_ms_level);
373
374
375 // Arrival time at half the range.
376
377 qualified_mass_spectrum.setDtInMilliSeconds(
378 tims_frame_csp.get()->getDriftTimeInMilliseconds(
379 mobility_scan_num_range_begin + (mobility_scan_num_range_width / 2)));
380
381 // 1/K0
382 qDebug() << "mobility_one_over_k0:" << mobility_one_over_k0
383 << "mobility_one_over_k0_range_begin:"
384 << mobility_one_over_k0_range_begin
385 << "mobility_one_over_k0_range_end"
386 << mobility_one_over_k0_range_end;
387
388 if(mobility_one_over_k0 == std::numeric_limits<double>::max() ||
389 mobility_one_over_k0_range_begin ==
390 std::numeric_limits<double>::max() ||
391 mobility_one_over_k0_range_end == std::numeric_limits<double>::max())
392 throw(
393 ExceptionNotPossible("Not possible that mobility_one_over_k0 and its "
394 "range are undefined."));
395
396 qualified_mass_spectrum.setParameterValue(
398 qualified_mass_spectrum.setParameterValue(
400 mobility_one_over_k0_range_begin);
401 qualified_mass_spectrum.setParameterValue(
403 mobility_one_over_k0_range_end);
404
405 // qDebug() << "mobility_scan_num_range_begin:"
406 // << mobility_scan_num_range_begin
407 // << "mobility_scan_num_range_end:" <<
408 // mobility_scan_num_range_end;
409
410 if(mobility_scan_num_range_begin == std::numeric_limits<quint32>::max() ||
411 mobility_scan_num_range_end == std::numeric_limits<quint32>::max())
413 "Not possible that mobility_scan_num_range values are undefined."));
414
415 qualified_mass_spectrum.setParameterValue(
417 mobility_scan_num_range_begin + (mobility_scan_num_range_width / 2));
418 qualified_mass_spectrum.setParameterValue(
420 mobility_scan_num_range_begin);
421 qualified_mass_spectrum.setParameterValue(
423 mobility_scan_num_range_end);
424
425 qualified_mass_spectrum.setParameterValue(
427 static_cast<qlonglong>(tims_frame_csp->getTotalNumberOfScans()));
428
429
430 Trace trace;
431
432 if(config.needPeakList())
433 {
434 // Provide these two variables for the function below to fill in the
435 // values. We will need them later.
436 quint32 min_mz_index_out = 0;
437 quint32 max_mz_index_out = 0;
438
439 if(asked_mz_range)
440 {
441 trace =
443 mz_index_merge_window,
444 mz_range_begin,
445 mz_range_end,
446 mobility_scan_num_range_begin,
447 mobility_scan_num_range_end,
448 min_mz_index_out,
449 max_mz_index_out);
450 }
451 else
452 {
453 trace =
455 mz_index_merge_window,
456 mobility_scan_num_range_begin,
457 mobility_scan_num_range_end,
458 min_mz_index_out,
459 max_mz_index_out);
460 }
461
462 // qDebug() << "Got min_mz_index_out:" << min_mz_index_out;
463 // qDebug() << "Got max_mz_index_out:" << max_mz_index_out;
464
465 qualified_mass_spectrum.setParameterValue(
467 min_mz_index_out);
468 qualified_mass_spectrum.setParameterValue(
470 max_mz_index_out);
471
472 qualified_mass_spectrum.setMassSpectrumSPtr(
473 std::make_shared<MassSpectrum>(trace));
474 qualified_mass_spectrum.setEmptyMassSpectrum(false);
475 }
476 else
477 {
478 qualified_mass_spectrum.setEmptyMassSpectrum(true);
479 }
480
481 handler.setQualifiedMassSpectrum(qualified_mass_spectrum);
482 }
483}
484
485
486void
488 [[maybe_unused]] SpectrumCollectionHandlerInterface &handler,
489 [[maybe_unused]] unsigned int ms_level)
490{
491 qDebug();
492}
493
494
495std::size_t
500
501
502
503Trace
505{
506
507 // We want to compute the TIC chromatogram, not load the chromatogram that
508 // is located in the SQL database.
509 //
510 // For this, we need to iterated into the frames and ask for MS1 spectra
511 // only. msp_timsData has that information:
512 //
513 // std::vector<FrameIdDescr> m_frameIdDescrList;
514 //
515 // and
516
517 // struct FrameIdDescr
518 // {
519 // std::size_t m_frameId; // frame id
520 // std::size_t m_size; // frame size (number of TOF scans in frame)
521 // std::size_t m_cumulSize; // cumulative size
522 // };
523
524 Trace tic_chromatogram;
525
526 const std::vector<FrameIdDescr> frame_descr_list =
528
529 for(FrameIdDescr frame_id_descr : frame_descr_list)
530 {
531 TimsFrameCstSPtr tims_frame_csp =
532 msp_timsData->getTimsFrameCstSPtrCached(frame_id_descr.m_frameId);
533 std::size_t scan_begin = 0;
534 std::size_t scan_end = tims_frame_csp->getTotalNumberOfScans() - 1;
535
536 // By convention, a TIC chromatogram is only performed using MS1
537 // spectra.
538 if(tims_frame_csp->getMsLevel() == 1)
539 {
540
541 // Retention times are in seconds in the Bruker world.
542 double rt = tims_frame_csp->getRtInSeconds();
543
544 tic_chromatogram.append(
546 tims_frame_csp->cumulateScanRangeIntensities(scan_begin,
547 scan_end)));
548 }
549 else
550 continue;
551 }
552
553 return tic_chromatogram;
554}
void setNativeId(const QString &native_id)
void setMsRunId(MsRunIdCstSPtr other)
void setSpectrumIndex(std::size_t index)
const QVariant getParameterValue(MsRunReadConfigParameter parameter) const
bool acceptMsLevel(std::size_t ms_level) const
bool acceptRetentionTimeInSeconds(double retention_time_in_seconds) const
const MsRunIdCstSPtr & getMsRunId() const
Class representing a fully specified mass spectrum.
void setDtInMilliSeconds(pappso_double rt)
Set the drift time in milliseconds.
void setMassSpectrumId(const MassSpectrumId &iD)
Set the MassSpectrumId.
void setMsLevel(uint ms_level)
Set the mass spectrum level.
void setParameterValue(QualifiedMassSpectrumParameter parameter, const QVariant &value)
void setMassSpectrumSPtr(MassSpectrumSPtr massSpectrum)
Set the MassSpectrumSPtr.
void setRtInSeconds(pappso_double rt)
Set the retention time in seconds.
void setEmptyMassSpectrum(bool is_empty_mass_spectrum)
interface to collect spectrums from the MsRunReader class
virtual void setQualifiedMassSpectrum(const QualifiedMassSpectrum &spectrum)=0
std::size_t getTotalScanCount() const
Definition timsdata.cpp:700
const std::vector< TimsFrameRecord > & getTimsFrameRecordList() const
TimsFrameCstSPtr getTimsFrameCstSPtrCached(std::size_t timsId)
get a Tims frame with his database ID but look in the cache first
Definition timsdata.cpp:953
const std::vector< FrameIdDescr > & getFrameIdDescrList() const
void getQualifiedMassSpectrumByGlobalScanIndex(const MsRunIdCstSPtr &msrun_id, QualifiedMassSpectrum &mass_spectrum, std::size_t global_scan_index, bool want_binary_data)
Definition timsdata.cpp:734
pappso::MassSpectrumCstSPtr getMassSpectrumCstSPtrByGlobalScanIndex(std::size_t index)
Definition timsdata.cpp:410
virtual std::size_t getTotalNumberOfScans() const
get the number of scans contained in this frame each scan represents an ion mobility slice
double getDriftTimeInMilliseconds(std::size_t scan_index) const
get drift time of a scan number in milliseconds
std::size_t getScanIndexFromOneOverK0(double one_over_k0) const
get the scan number from a given 1/Ko mobility value
double getRtInSeconds() const
unsigned int getMsLevel() const
double getOneOverK0Transformation(std::size_t scan_index) const
get 1/K0 value of a given scan (mobility value)
virtual Trace combineScansToTraceWithDowngradedMzResolution(std::size_t mzindex_merge_window, std::size_t scanNumBegin, std::size_t scanNumEnd, quint32 &mz_minimum_index, quint32 &mz_maximum_index) const override
cumulate spectrum given a scan number range need the binary file The intensities are normalized with ...
virtual quint64 cumulateScanRangeIntensities(std::size_t scanNumBegin, std::size_t scanNumEnd) const override
...
virtual Trace combineScansToTraceWithDowngradedMzResolution2(std::size_t mz_index_merge_window, double mz_range_begin, double mz_range_end, std::size_t mobility_scan_begin, std::size_t mobility_scan_end, quint32 &mz_minimum_index_out, quint32 &mz_maximum_index_out) const override
cumulate spectrum given a scan number range need the binary file The intensities are normalized with ...
virtual std::size_t spectrumListSize() const override
get the totat number of spectrum conained in the MSrun data file
virtual void readSpectrumCollection(SpectrumCollectionHandlerInterface &handler) override
function to visit an MsRunReader and get each Spectrum in a spectrum collection handler
virtual MassSpectrumCstSPtr massSpectrumCstSPtr(std::size_t spectrum_index) override
virtual QualifiedMassSpectrum qualifiedMassSpectrum(std::size_t spectrum_index, bool want_binary_data=true) const override
get a QualifiedMassSpectrum class given its scan number
virtual MassSpectrumSPtr massSpectrumSPtr(std::size_t spectrum_index) override
get a MassSpectrumSPtr class given its spectrum index
virtual void readSpectrumCollection2(const MsRunReadConfig &config, SpectrumCollectionHandlerInterface &handler) override
TimsFramesMsRunReader(MsRunIdCstSPtr &msrun_id_csp)
virtual void readSpectrumCollectionByMsLevel(SpectrumCollectionHandlerInterface &handler, unsigned int ms_level) override
function to visit an MsRunReader and get each Spectrum in a spectrum collection handler by Ms Levels
A simple container of DataPoint instances.
Definition trace.h:148
size_t append(const DataPoint &data_point)
appends a datapoint and return new size
Definition trace.cpp:648
process interrupted exception
tries to keep as much as possible monoisotopes, removing any possible C13 peaks and changes multichar...
Definition aa.cpp:39
std::shared_ptr< const MsRunId > MsRunIdCstSPtr
Definition msrunid.h:46
std::shared_ptr< const MassSpectrum > MassSpectrumCstSPtr
@ TimsFrameMzIndexBegin
Bruker's timsTOF mz index frame start range.
@ IonMobOneOverK0Begin
1/K0 range's begin value
@ IonMobOneOverK0End
1/K0 range's end value
@ TimsFrameScansCount
Bruker's timsTOF frame's total ion mobility slots.
@ TimsFrameMzIndexEnd
Bruker's timsTOF mz index frame end range.
@ rt
Retention time.
std::shared_ptr< MassSpectrum > MassSpectrumSPtr
std::shared_ptr< const TimsFrame > TimsFrameCstSPtr
Definition timsframe.h:43
std::size_t m_globalScanIndex
Definition timsdata.h:63
handle specific data for DDA MS runs