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
29#include "../../exception/exceptionnotimplemented.h"
30#include "../../exception/exceptioninterrupted.h"
31#include "../../exception/exceptionnotpossible.h"
32#include <QDebug>
33
34using namespace pappso;
35
37 : MsRunReader(msrun_id_csp)
38{
39 qDebug() << "Now initializing the TimsFramesMsRunReader.";
40
41 initialize();
42}
43
48
49void
51{
52 msp_timsData = std::make_shared<TimsData>(mcsp_msRunId.get()->getFileName());
53 if(msp_timsData == nullptr)
54 {
55 throw PappsoException(
56 QObject::tr("ERROR in TimsFramesMsRunReader::initialize "
57 "msp_timsData is null for MsRunId %1")
58 .arg(mcsp_msRunId.get()->toString()));
59 }
60}
61
62
63bool
64TimsFramesMsRunReader::accept(const QString &file_name) const
65{
66 qDebug() << file_name;
67 return true;
68}
69
70
73 [[maybe_unused]] std::size_t spectrum_index)
74{
76 QObject::tr("Not yet implemented in TimsFramesMsRunReader %1.\n")
77 .arg(__LINE__));
78
80}
81
82
85{
86 return msp_timsData->getMassSpectrumCstSPtrByRawIndex(spectrum_index);
87}
88
89
92 bool want_binary_data) const
93{
94
95 QualifiedMassSpectrum mass_spectrum;
96
97 msp_timsData->getQualifiedMassSpectrumByRawIndex(
98 getMsRunId(), mass_spectrum, spectrum_index, want_binary_data);
99 return mass_spectrum;
100}
101
102
103void
106{
107 // qDebug() << "Reading the spectrum collection with no specific
108 // configuration.";
109 MsRunReadConfig config;
110 readSpectrumCollection2(config, handler);
111}
112
113
114void
117{
118 // qDebug().noquote() << "Reading the spectrum collection with this "
119 // "specific configuration:"
120 // << config.toString();
121
122 std::vector<std::size_t> subset_of_tims_frame_ids;
123
124
125 bool asked_ion_mobility_scan_num_range = false;
126
127 quint32 mobility_scan_num_range_begin = std::numeric_limits<quint32>::max();
128 quint32 mobility_scan_num_range_end = std::numeric_limits<quint32>::max();
129 quint32 mobility_scan_num_range_width = std::numeric_limits<quint32>::max();
130
131 double mobility_one_over_k0 = std::numeric_limits<double>::max();
132 double mobility_one_over_k0_range_begin = std::numeric_limits<double>::max();
133 double mobility_one_over_k0_range_end = std::numeric_limits<double>::max();
134
135 if(!config
136 .getParameterValue(
138 .isNull() &&
139 !config
140 .getParameterValue(
142 .isNull())
143 {
144 mobility_scan_num_range_begin =
145 config
148 .toUInt();
149 mobility_scan_num_range_end =
150 config
153 .toUInt();
154
155 // We need the range width below.
156 mobility_scan_num_range_width =
157 mobility_scan_num_range_end + 1 - mobility_scan_num_range_begin;
158
159 asked_ion_mobility_scan_num_range = true;
160
161 // Be sure to check in the frames loop below that the user might
162 // have asked for an ion mobility range but on the basis of the 1/K0 unit.
163 }
164
165 const std::vector<FrameIdDescr> &frame_id_descr_list =
166 msp_timsData->getFrameIdDescrList();
167
168 // Just for the feedback to the user.
169 std::size_t scan_count = 0;
170
171 for(auto const &frame_record : msp_timsData->getTimsFrameRecordList())
172 {
173 if(handler.shouldStop())
174 {
175 // qDebug() << "The operation was cancelled. Breaking the loop.";
177 QObject::tr("Reading timsTOF data cancelled by the user."));
178 }
179
180 if(frame_record.frame_id == 0)
181 continue;
182
183 if(!config.acceptRetentionTimeInSeconds(frame_record.frame_time))
184 continue;
185
186 std::size_t ms_level = 2;
187 if(frame_record.msms_type == 0)
188 ms_level = 1;
189
190 if(!config.acceptMsLevel(ms_level))
191 continue;
192
193 subset_of_tims_frame_ids.push_back(frame_record.frame_id);
194
195 if(mobility_scan_num_range_width != std::numeric_limits<int>::max())
196 {
197 scan_count += mobility_scan_num_range_width;
198 }
199 else
200 {
201 scan_count += frame_id_descr_list[frame_record.frame_id].m_size;
202 }
203 }
204
205 // At this point, we have a subset of frame records.
206 std::size_t frame_count = subset_of_tims_frame_ids.size();
207 qDebug() << "The number of retained RT range- and MS level-matching frames : "
208 << frame_count;
209
210 // Inform the handler of the spectrum list so that it can handle feedback to
211 // the user.
212
213 // FIXME:
214 // Either we document the number of frames (because we assume we will
215 // flatten them all)...
216 handler.spectrumListHasSize(frame_count);
217 // Or we document the number of actual scans because we might not flatten
218 // all the frames.
219 // handler.spectrumListHasSize(scan_count);
220
221 // Check for m/z range selection
222 bool asked_mz_range = false;
223 double mz_range_begin = -1;
224 double mz_range_end = -1;
225
227 {
228 asked_mz_range = true;
229
230 mz_range_begin =
232 .toDouble();
233 mz_range_end =
235 .toDouble();
236
237 // qDebug() << "The m/z range asked is: " << mz_range_begin << "--"
238 // << mz_range_end;
239 }
240
241 // Check for m/z resolution downgrading
242 // The idea is that we merge a number of mz indices into a single index,
243 // which is essentially an increase of the m/z bin size, and therefore
244 // of the resolution/definition of the mass spectrum.
245 std::size_t mz_index_merge_window = 0;
246 if(!config
247 .getParameterValue(
249 .isNull())
250 {
251 mz_index_merge_window =
252 config
255 .toUInt();
256
257 // qDebug() << "mz_index_merge_window=" << mz_index_merge_window;
258 }
259
260 // The scan index is the index of the scan in the *whole* mass data file, it
261 // is a sequential number of scans over all the frames.
262 std::size_t scan_index = 0; // iterate in each spectrum
263
264 for(std::size_t tims_frame_id : subset_of_tims_frame_ids)
265 {
266 if(handler.shouldStop())
267 {
268 // qDebug() << "The operation was cancelled. Breaking the loop.";
270 QObject::tr("Reading timsTOF data cancelled by the user."));
271 }
272
273 // qDebug() << "tims_frame_id=" << tims_frame_id;
274
275 const FrameIdDescr &current_frame_record =
276 frame_id_descr_list[tims_frame_id];
277 // qDebug() << "tims_frame_id=" << tims_frame_id;
278
279 scan_index = current_frame_record.m_cumulSize;
280
281 TimsFrameCstSPtr tims_frame_csp =
282 msp_timsData->getTimsFrameCstSPtrCached(tims_frame_id);
283
284 // If the user wants to select 1/Ko values in a given range, we need to
285 // compute the ion mobility scan value starting from that 1/Ko value in
286 // *each* frame. Note that the computed mobility_scan_num_begin and
287 // mobility_scan_num_end would override thoses possibly set with
288 // TimsFramesMsRunReader_mobility_index_begin/end above.
289
290 if(!config
291 .getParameterValue(
293 .isNull() &&
294 !config
295 .getParameterValue(
297 .isNull())
298 {
299 mobility_one_over_k0_range_begin =
300 config
303 .toDouble();
304
305 mobility_one_over_k0_range_end =
306 config
309 .toDouble();
310
311 mobility_scan_num_range_begin =
312 tims_frame_csp.get()->getScanNumFromOneOverK0(
313 mobility_one_over_k0_range_begin);
314
315 mobility_scan_num_range_end =
316 tims_frame_csp.get()->getScanNumFromOneOverK0(
317 mobility_one_over_k0_range_end);
318
319 asked_ion_mobility_scan_num_range = true;
320 }
321
322 // Now that we know if the user has asked for an ion mobility range,
323 // either using scan indices or 1/K0 values, we need to double check the
324 // range borders.
325 quint32 count_of_mobility_scans = tims_frame_csp->getTotalNumberOfScans();
326
327 if(asked_ion_mobility_scan_num_range)
328 {
329 if(mobility_scan_num_range_end > (count_of_mobility_scans - 1))
330 {
331 mobility_scan_num_range_end = count_of_mobility_scans - 1;
332 }
333 }
334 else
335 {
336 mobility_scan_num_range_begin = 0;
337 mobility_scan_num_range_end = count_of_mobility_scans - 1;
338 }
339
340 // Now that we know the mobility index range, if we did not set the
341 // mobility one over K0 because that was not the unit used by
342 // the caller, then we can compute these values and set them
343 // later to the qualified mass spectrum parameters.
344 if(mobility_one_over_k0_range_begin == std::numeric_limits<double>::max())
345 mobility_one_over_k0_range_begin =
346 tims_frame_csp->getOneOverK0Transformation(
347 mobility_scan_num_range_begin);
348 if(mobility_one_over_k0_range_end == std::numeric_limits<double>::max())
349 mobility_one_over_k0_range_end =
350 tims_frame_csp->getOneOverK0Transformation(
351 mobility_scan_num_range_end);
352
353 mobility_scan_num_range_width =
354 mobility_scan_num_range_end + 1 - mobility_scan_num_range_begin;
355
356 // We want to provide the inverse mobility for the scan that sits in the
357 // middle of the defined range or the whole range if none is defined..
358 mobility_one_over_k0 = tims_frame_csp.get()->getScanNumFromOneOverK0(
359 mobility_scan_num_range_begin + (mobility_scan_num_range_width / 2));
360
361 // Now, with or without the peak list, we have to craft a qualified mass
362 // spectrum that will hold all the data about the data in it.
363 QualifiedMassSpectrum mass_spectrum;
364
365 MassSpectrumId spectrum_id;
366
367 spectrum_id.setSpectrumIndex(tims_frame_id);
368 spectrum_id.setMsRunId(getMsRunId());
369
370 // Can be modified to add bits that might help our case
371 spectrum_id.setNativeId(
372 QString("frame_id=%1 global_scan_index=%2 im_scan_range_begin=%3 "
373 "im_scan_range_end=%4")
374 .arg(tims_frame_id)
375 .arg(scan_index)
376 .arg(mobility_scan_num_range_begin)
377 .arg(mobility_scan_num_range_end));
378
379 mass_spectrum.setMassSpectrumId(spectrum_id);
380
381 // We want to document the retention time!
382 mass_spectrum.setRtInSeconds(tims_frame_csp.get()->getTime());
383
384 // We do want to document the ms level of the spectrum and possibly
385 // the precursor's m/z and charge.
386 unsigned int frame_ms_level = tims_frame_csp.get()->getMsLevel();
387 mass_spectrum.setMsLevel(frame_ms_level);
388
389 // Arrival time at half the range.
390 mass_spectrum.setDtInMilliSeconds(tims_frame_csp.get()->getDriftTime(
391 mobility_scan_num_range_begin + (mobility_scan_num_range_width / 2)));
392
393 // 1/K0
394 qDebug() << "mobility_one_over_k0:" << mobility_one_over_k0
395 << "mobility_one_over_k0_range_begin:"
396 << mobility_one_over_k0_range_begin
397 << "mobility_one_over_k0_range_end"
398 << mobility_one_over_k0_range_end;
399
400 if(mobility_one_over_k0 == std::numeric_limits<double>::max() ||
401 mobility_one_over_k0_range_begin ==
402 std::numeric_limits<double>::max() ||
403 mobility_one_over_k0_range_end == std::numeric_limits<double>::max())
404 throw(
405 ExceptionNotPossible("Not possible that mobility_one_over_k0 and its "
406 "range are undefined."));
407
408 mass_spectrum.setParameterValue(
410 mobility_one_over_k0);
411 mass_spectrum.setParameterValue(
413 mobility_one_over_k0_range_begin);
414 mass_spectrum.setParameterValue(
416 mobility_one_over_k0_range_end);
417
418 // qDebug() << "mobility_scan_num_range_begin:"
419 // << mobility_scan_num_range_begin
420 // << "mobility_scan_num_range_end:" <<
421 // mobility_scan_num_range_end;
422
423 if(mobility_scan_num_range_begin == std::numeric_limits<quint32>::max() ||
424 mobility_scan_num_range_end == std::numeric_limits<quint32>::max())
426 "Not possible that mobility_scan_num_range values are undefined."));
427
428 mass_spectrum.setParameterValue(
430 mobility_scan_num_range_begin + (mobility_scan_num_range_width / 2));
431 mass_spectrum.setParameterValue(
433 mobility_scan_num_range_begin);
434 mass_spectrum.setParameterValue(
436 mobility_scan_num_range_end);
437
438 mass_spectrum.setParameterValue(
440 static_cast<qlonglong>(tims_frame_csp->getTotalNumberOfScans()));
441
442 Trace trace;
443
444 if(config.needPeakList())
445 {
446 // Provide these two variables for the function below to fill in the
447 // values. We will need them later.
448 quint32 min_mz_index_out = 0;
449 quint32 max_mz_index_out = 0;
450
451 if(asked_mz_range)
452 {
453 trace = tims_frame_csp->cumulateScansToTraceMzDownResolution2(
454 mz_index_merge_window,
455 mz_range_begin,
456 mz_range_end,
457 mobility_scan_num_range_begin,
458 mobility_scan_num_range_end,
459 min_mz_index_out,
460 max_mz_index_out);
461 }
462 else
463 {
464 trace = tims_frame_csp->cumulateScansToTraceMzDownResolution(
465 mz_index_merge_window,
466 mobility_scan_num_range_begin,
467 mobility_scan_num_range_end,
468 min_mz_index_out,
469 max_mz_index_out);
470 }
471
472 // qDebug() << "Got min_mz_index_out:" << min_mz_index_out;
473 // qDebug() << "Got max_mz_index_out:" << max_mz_index_out;
474
475 mass_spectrum.setParameterValue(
477 min_mz_index_out);
478 mass_spectrum.setParameterValue(
480 max_mz_index_out);
481
482 mass_spectrum.setMassSpectrumSPtr(
483 std::make_shared<MassSpectrum>(trace));
484 mass_spectrum.setEmptyMassSpectrum(false);
485 }
486 else
487 {
488 mass_spectrum.setEmptyMassSpectrum(true);
489 }
490
491 handler.setQualifiedMassSpectrum(mass_spectrum);
492 }
493}
494
495
496void
498 [[maybe_unused]] SpectrumCollectionHandlerInterface &handler,
499 [[maybe_unused]] unsigned int ms_level)
500{
501 qDebug();
502}
503
504
505std::size_t
507{
508 return msp_timsData->getTotalNumberOfScans();
509}
510
511
512bool
514{
515 return false;
516}
517
518
519bool
521{
522 msp_timsData = nullptr;
523 return true;
524}
525
526bool
528{
529 if(msp_timsData == nullptr)
530 {
531 initialize();
532 }
533 return true;
534}
535
536
539 std::size_t spectrum_index [[maybe_unused]],
540 pappso::PrecisionPtr precision [[maybe_unused]]) const
541{
542 throw ExceptionNotImplemented(QObject::tr("Not implemented %1 %2 %3")
543 .arg(__FILE__)
544 .arg(__FUNCTION__)
545 .arg(__LINE__));
546}
547
550 const pappso::QualifiedMassSpectrum &mass_spectrum [[maybe_unused]],
551 pappso::PrecisionPtr precision [[maybe_unused]]) const
552{
553 throw ExceptionNotImplemented(QObject::tr("Not implemented %1 %2 %3")
554 .arg(__FILE__)
555 .arg(__FUNCTION__)
556 .arg(__LINE__));
557}
558
565
566
567Trace
569{
570 // Use the Sqlite database to fetch the total ion current chromatogram (TIC
571 // chromatogram).
572
574
575 return msp_timsData->getTicChromatogram();
576}
577
578
579Trace
581{
582
583 // We want to compute the TIC chromatogram, not load the chromatogram that
584 // is located in the SQL database.
585 //
586 // For this, we need to iterated into the frames and ask for MS1 spectra
587 // only. msp_timsData has that information:
588 //
589 // std::vector<FrameIdDescr> m_frameIdDescrList;
590 //
591 // and
592
593 // struct FrameIdDescr
594 // {
595 // std::size_t m_frameId; // frame id
596 // std::size_t m_size; // frame size (number of TOF scans in frame)
597 // std::size_t m_cumulSize; // cumulative size
598 // };
599
600 Trace tic_chromatogram;
601
602 const std::vector<FrameIdDescr> frame_descr_list =
603 msp_timsData->getFrameIdDescrList();
604
605 for(FrameIdDescr frame_id_descr : frame_descr_list)
606 {
607 TimsFrameCstSPtr tims_frame_csp =
608 msp_timsData->getTimsFrameCstSPtrCached(frame_id_descr.m_frameId);
609 std::size_t scan_begin = 0;
610 std::size_t scan_end = tims_frame_csp->getTotalNumberOfScans() - 1;
611
612 // By convention, a TIC chromatogram is only performed using MS1
613 // spectra.
614 if(tims_frame_csp->getMsLevel() == 1)
615 {
616
617 // Retention times are in seconds in the Bruker world.
618 double rt = tims_frame_csp->getTime();
619
620 tic_chromatogram.append(DataPoint(
621 rt,
622 tims_frame_csp->cumulateScansIntensities(scan_begin, scan_end)));
623 }
624 else
625 continue;
626 }
627
628 return tic_chromatogram;
629}
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
base class to read MSrun the only way to build a MsRunReader object is to use the MsRunReaderFactory
Definition msrunreader.h:63
MsRunIdCstSPtr mcsp_msRunId
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
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 pappso::XicCoordSPtr newXicCoordSPtrFromQualifiedMassSpectrum(const pappso::QualifiedMassSpectrum &mass_spectrum, pappso::PrecisionPtr precision) const override
get a xic coordinate object from a given spectrum
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 Trace getTicChromatogram() override
get a TIC chromatogram
virtual bool releaseDevice() override
release data back end device if a the data back end is released, the developper has to use acquireDev...
virtual void readSpectrumCollection2(const MsRunReadConfig &config, SpectrumCollectionHandlerInterface &handler) override
virtual TimsDataSp getTimsDataSPtr()
give an access to the underlying raw data pointer
virtual bool hasScanNumbers() const override
tells if spectra can be accessed using scan numbers by default, it returns false. Only overrided func...
virtual bool acquireDevice() override
acquire data back end device
virtual bool accept(const QString &file_name) const override
tells if the reader is able to handle this file must be implemented by private MS run reader,...
TimsFramesMsRunReader(MsRunIdCstSPtr &msrun_id_csp)
virtual pappso::XicCoordSPtr newXicCoordSPtrFromSpectrumIndex(std::size_t spectrum_index, pappso::PrecisionPtr precision) const override
get a xic coordinate object from a given spectrum index
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
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< TimsData > TimsDataSp
shared pointer on a TimsData object
Definition timsdata.h:50
std::shared_ptr< const MassSpectrum > MassSpectrumCstSPtr
@ TimsIonMobScanOneOverK0
1/kO of a simple scan
@ TimsFrameMzIndexBegin
Bruker's timsTOF mz index frame start range.
@ 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::shared_ptr< XicCoord > XicCoordSPtr
Definition xiccoord.h:43
std::size_t m_cumulSize
Definition timsdata.h:57