libpappsomspp
Library for mass spectrometry
Loading...
Searching...
No Matches
timsmsrunreader.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
28#include "timsmsrunreader.h"
33#include <QDebug>
34
35using namespace pappso;
36
38 : TimsMsRunReaderBase(msrun_id_csp)
39{
40 initialize();
41}
43 : TimsMsRunReaderBase(msrun_reader_base)
44{
45 initialize();
46}
47
48
52
53
55TimsMsRunReader::massSpectrumSPtr([[maybe_unused]] std::size_t spectrum_index)
56{
58 QObject::tr("Not yet implemented in TimsMsRunReader %1.\n").arg(__LINE__));
60}
61
62
64TimsMsRunReader::massSpectrumCstSPtr(std::size_t spectrum_index)
65{
67}
68
69
71TimsMsRunReader::qualifiedMassSpectrum(std::size_t spectrum_index,
72 bool want_binary_data) const
73{
74
75 QualifiedMassSpectrum mass_spectrum;
76
78 getMsRunId(), mass_spectrum, spectrum_index, want_binary_data);
79 return mass_spectrum;
80}
81
82
83void
89
90void
93{
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 bool asked_ion_mobility_scan_num_range = false;
102
103 quint32 mobility_scan_num_range_begin =
104 std::numeric_limits<quint32>::quiet_NaN();
105 quint32 mobility_scan_num_range_end =
106 std::numeric_limits<quint32>::quiet_NaN();
107 quint32 mobility_scan_num_range_width =
108 std::numeric_limits<quint32>::quiet_NaN();
109
110 double mobility_one_over_k0_range_begin =
111 std::numeric_limits<double>::quiet_NaN();
112 double mobility_one_over_k0_range_end =
113 std::numeric_limits<double>::quiet_NaN();
114
115 if(!config
116 .getParameterValue(
118 .isNull() &&
119 !config
120 .getParameterValue(
122 .isNull())
123 {
124 mobility_scan_num_range_begin =
125 config
128 .toUInt();
129 mobility_scan_num_range_end =
130 config
133 .toUInt();
134
135 // We need the range width below.
136 mobility_scan_num_range_width =
137 mobility_scan_num_range_end + 1 - mobility_scan_num_range_begin;
138
139 asked_ion_mobility_scan_num_range = true;
140
141 // Be sure to check in the frames loop below that the user might
142 // have asked for an ion mobility range but on the basis of the 1/K0 unit.
143 }
144
145 const std::vector<FrameIdDescr> &frame_id_descr_list =
147
148 // Just for the feedback to the user.
149 std::size_t scan_count = 0;
150
151 for(auto const &frame_record : msp_timsData->getTimsFrameRecordList())
152 {
153 if(handler.shouldStop())
154 {
155 // qDebug() << "The operation was cancelled. Breaking the loop.";
157 QObject::tr("Reading timsTOF data cancelled by the user."));
158 }
159
160 if(frame_record.frame_id == 0)
161 continue;
162
163 if(!config.acceptRetentionTimeInSeconds(frame_record.frame_time))
164 continue;
165
166 std::size_t ms_level = 2;
167 if(frame_record.msms_type == 0)
168 ms_level = 1;
169
170 if(!config.acceptMsLevel(ms_level))
171 continue;
172
173 subset_of_tims_frame_ids.push_back(frame_record.frame_id);
174
175 if(mobility_scan_num_range_width)
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 // Inform the handler of the spectrum list so that it can handle feedback to
190 // the user.
191 handler.spectrumListHasSize(scan_count);
192
193 // Check for m/z range selection
194 double mz_range_begin = -1;
195 double mz_range_end = -1;
196
198 .isNull() &&
200 {
201 mz_range_begin =
203 .toDouble();
204
205 mz_range_end =
207 .toDouble();
208
209 // qDebug() << "The m/z range asked is: " << mz_range_begin
210 // << "--" << mz_range_end;
211 }
212
213 // Check for m/z resolution downgrading (mz bins merge)
214 // The idea is that we merge a number of mz indices into a single index,
215 // which is essentially an increase of the m/z bin size, and therefore
216 // a reduction of the resolution/definition of the mass spectrum.
217 std::size_t mz_index_merge_window = 0;
218 if(!config
219 .getParameterValue(
221 .isNull())
222 {
223 mz_index_merge_window =
224 config
227 .toUInt();
228
229 // qDebug() << "mz_index_merge_window=" << mz_index_merge_window;
230 }
231
232 std::size_t number_of_mobility_scans_set_as_qualified_mass_spectra = 0;
233
234 for(std::size_t tims_frame_id : subset_of_tims_frame_ids)
235 {
236 if(handler.shouldStop())
237 {
238 // qDebug() << "The operation was cancelled. Breaking the loop.";
240 QObject::tr("Reading timsTOF data cancelled by the user."));
241 }
242
243 qDebug() << "tims_frame_id=" << tims_frame_id;
244
245 const FrameIdDescr &current_frame_record =
246 frame_id_descr_list[tims_frame_id];
247
248 TimsFrameCstSPtr tims_frame_csp =
250
251 qDebug() << "tims_frame_id=" << tims_frame_id;
252
253 // If the user wants to select 1/Ko values in a given range, we need to
254 // compute the ion mobility scan value starting from that 1/Ko value in
255 // *each* frame. Note that the computed mobility_scan_num_begin and
256 // mobility_scan_num_end would override thoses possibly set with
257 // TimsFramesMsRunReader_mobility_index_begin/end above.
258
259 if(!config
260 .getParameterValue(
262 .isNull() &&
263 !config
264 .getParameterValue(
266 .isNull())
267 {
268 mobility_one_over_k0_range_begin =
269 config
272 .toDouble();
273
274 mobility_one_over_k0_range_end =
275 config
278 .toDouble();
279
280 mobility_scan_num_range_begin =
281 tims_frame_csp.get()->getScanIndexFromOneOverK0(
282 mobility_one_over_k0_range_begin);
283
284 mobility_scan_num_range_end =
285 tims_frame_csp.get()->getScanIndexFromOneOverK0(
286 mobility_one_over_k0_range_end);
287
288 asked_ion_mobility_scan_num_range = true;
289 }
290
291 qDebug() << "tims_frame_id=" << tims_frame_id;
292 // Now that we know if the user has asked for an ion mobility range,
293 // either using scan indices or 1/K0 values, we need to double check the
294 // range borders.
295
296 quint32 count_of_mobility_scans = tims_frame_csp->getTotalNumberOfScans();
297
298 if(asked_ion_mobility_scan_num_range)
299 {
300 if(mobility_scan_num_range_end > (count_of_mobility_scans - 1))
301 {
302 mobility_scan_num_range_end = count_of_mobility_scans - 1;
303 }
304 }
305 else
306 {
307 mobility_scan_num_range_begin = 0;
308 mobility_scan_num_range_end = count_of_mobility_scans - 1;
309 }
310
311 qDebug() << "tims_frame_id=" << tims_frame_id;
312 // Now, with or without the peak list, we have to craft a qualified mass
313 // spectrum that will hold all the data about the data in it.
314 QualifiedMassSpectrum mass_spectrum;
315
316 MassSpectrumId spectrum_id;
317
318 spectrum_id.setSpectrumIndex(tims_frame_id);
319 spectrum_id.setMsRunId(getMsRunId());
320
321 mass_spectrum.setMassSpectrumId(spectrum_id);
322
323 qDebug() << "tims_frame_id=" << tims_frame_id;
324 // We want to document the retention time!
325 mass_spectrum.setRtInSeconds(tims_frame_csp.get()->getRtInSeconds());
326
327 // We do want to document the ms level of the spectrum and possibly
328 // the precursor's m/z and charge.
329 unsigned int frame_ms_level = tims_frame_csp.get()->getMsLevel();
330 mass_spectrum.setMsLevel(frame_ms_level);
331
332 qDebug() << "tims_frame_id=" << tims_frame_id;
333 std::vector<TimsDdaPrecursors::SpectrumDescr> dda_spectrum_descr_list;
334 TimsDiaSlices::MsMsWindowGroup *p_dia_window_group = nullptr;
335 std::size_t frame_global_slice_begin = 0;
336 if(frame_ms_level > 1)
337 {
338 if(msp_timsData.get()->isDdaRun())
339 {
340 TimsDdaPrecursors *dda_precursors_p =
342 dda_spectrum_descr_list =
343 dda_precursors_p->getSpectrumDescrListByFrameId(tims_frame_id);
344 }
345
346 if(msp_timsData.get()->isDiaRun())
347 {
348
349 qDebug() << "tims_frame_id=" << tims_frame_id;
350 p_dia_window_group = msp_timsData.get()
353 .at(tims_frame_id);
354 frame_global_slice_begin =
355 msp_timsData.get()
357 ->getGlobalSliceIndexBeginByFrameId(tims_frame_id);
358 }
359 }
360
361 qDebug() << "tims_frame_id=" << tims_frame_id;
362 // The scan index is the index of the scan in the *whole* mass data file,
363 // it is a sequential number of scans over all the frames.
364 std::size_t scan_index = current_frame_record.m_globalScanIndex -
365 current_frame_record.m_scanCount +
366 mobility_scan_num_range_begin;
367
368 QualifiedMassSpectrum mass_spectrum_no_precursor_data(mass_spectrum);
369 for(quint32 iter_scan_index = mobility_scan_num_range_begin;
370 iter_scan_index <= mobility_scan_num_range_end;
371 iter_scan_index++)
372 {
373 mass_spectrum = mass_spectrum_no_precursor_data;
374 mass_spectrum.getMassSpectrumId().setSpectrumIndex(scan_index);
375
376 mass_spectrum.getMassSpectrumId().setNativeId(
377 QString("frame_id=%1 scan_index=%2 global_scan_index=%3")
378 .arg(tims_frame_id)
379 .arg(iter_scan_index)
380 .arg(scan_index));
381
382 // qDebug() << "iter:" << iter;
383
384 // qDebug() << "Scan's mz_minimum_index:" << mz_minimum_index_out
385 // << "and mz_maximum_index:" << mz_maximum_index_out;
386
387 // Arrival time
388 mass_spectrum.setDtInMilliSeconds(
389
390 tims_frame_csp.get()->getDriftTimeInMilliseconds(iter_scan_index));
391 // 1/K0
392 mass_spectrum.setParameterValue(
394 tims_frame_csp.get()->getOneOverK0Transformation(iter_scan_index));
395
396 if(dda_spectrum_descr_list.size() > 0)
397 {
398 std::size_t scan_index = (std::size_t)iter_scan_index;
399 auto it_spectrum_descr = std::find_if(
400 dda_spectrum_descr_list.begin(),
401 dda_spectrum_descr_list.end(),
402 [scan_index](
403 const TimsDdaPrecursors::SpectrumDescr &spectrum_descr) {
404 if(scan_index < spectrum_descr.scan_mobility_start)
405 return false;
406 if(scan_index > spectrum_descr.scan_mobility_end)
407 return false;
408 return true;
409 });
410
411 if(it_spectrum_descr != dda_spectrum_descr_list.end())
412 {
413 qDebug() << "scan_index=" << scan_index
414 << " spectrum_descr.scan_mobility_end="
415 << it_spectrum_descr->scan_mobility_end;
416 mass_spectrum.appendPrecursorIonData(
417 it_spectrum_descr->precursor_ion_data);
418
419 mass_spectrum.setPrecursorNativeId(
420 QString(
421 "frame_id=%1 begin=%2 end=%3 precursor=%4 idxms1=%5")
422 .arg(it_spectrum_descr->parent_frame)
423 .arg(it_spectrum_descr->scan_mobility_start)
424 .arg(it_spectrum_descr->scan_mobility_end)
425 .arg(it_spectrum_descr->precursor_id)
426 .arg(it_spectrum_descr->ms1_index));
427
428 mass_spectrum.setParameterValue(
430 it_spectrum_descr->isolationMz);
431 mass_spectrum.setParameterValue(
433 it_spectrum_descr->isolationWidth);
434
435 mass_spectrum.setParameterValue(
437 it_spectrum_descr->collisionEnergy);
438 mass_spectrum.setParameterValue(
440 (quint64)it_spectrum_descr->precursor_id);
441 }
442 }
443 if(p_dia_window_group != nullptr)
444 {
445
446
447 std::size_t scan_index = (std::size_t)iter_scan_index;
448 auto it_dia_window = std::find_if(
449 p_dia_window_group->begin(),
450 p_dia_window_group->end(),
451 [scan_index](const TimsDiaSlices::MsMsWindow &dia_window) {
452 if(scan_index < dia_window.ScanNumBegin)
453 return false;
454 if(scan_index > dia_window.ScanNumEnd)
455 return false;
456 return true;
457 });
458
459 if(it_dia_window != p_dia_window_group->end())
460 {
461 qDebug() << "scan_index=" << scan_index
462 << " it_dia_window->ScanNumEnd="
463 << it_dia_window->ScanNumEnd;
464
465 mass_spectrum.setPrecursorNativeId(
466 QString("window_group=%1 begin=%2 end=%3 frame=%4 scan=%5 "
467 "global_slice_id=%6")
468 .arg(it_dia_window->WindowGroup)
469 .arg(it_dia_window->ScanNumBegin)
470 .arg(it_dia_window->ScanNumEnd)
471 .arg(msp_timsData.get()
473 ->getLastMs1FrameIdByMs2FrameId(tims_frame_id))
474 .arg(iter_scan_index)
475 .arg(it_dia_window->SliceIndex +
476 frame_global_slice_begin));
477
478 mass_spectrum.setParameterValue(
480 it_dia_window->IsolationMz);
481 mass_spectrum.setParameterValue(
483 it_dia_window->IsolationWidth);
484
485 mass_spectrum.setParameterValue(
487 it_dia_window->CollisionEnergy);
488 }
489 }
490 if(config.needPeakList())
491 {
492 quint32 mz_minimum_index_out = 0;
493 quint32 mz_maximum_index_out = 0;
494
495 auto raw_trace =
496 tims_frame_csp.get()->getMobilityScan(iter_scan_index,
497 mz_index_merge_window,
498 mz_range_begin,
499 mz_range_end,
500 mz_minimum_index_out,
501 mz_maximum_index_out);
502
503 // qDebug() << "Ion mobility scan's raw trace size:"
504 // << raw_trace.size();
505
506 mass_spectrum.setEmptyMassSpectrum(false);
507
508 mass_spectrum.setParameterValue(
510 mz_minimum_index_out);
511 mass_spectrum.setParameterValue(
513 mz_maximum_index_out);
514
515
516 qDebug();
517 mass_spectrum.setMassSpectrumSPtr(
518 std::make_shared<MassSpectrum>(raw_trace));
519
520 qDebug() << mass_spectrum.getRtInSeconds();
521 }
522 else
523 {
524 mass_spectrum.setEmptyMassSpectrum(true);
525 }
526
527 // qDebug() << "mz_index_merge_window=" << mz_index_merge_window;
528 handler.setQualifiedMassSpectrum(mass_spectrum);
529 ++number_of_mobility_scans_set_as_qualified_mass_spectra;
530 scan_index++;
531 }
532 }
533
534 qDebug() << "Total number of loaded mass spectra:"
535 << number_of_mobility_scans_set_as_qualified_mass_spectra;
536}
537
538void
540 SpectrumCollectionHandlerInterface &handler, unsigned int ms_level)
541{
542
543 qDebug();
544
545 try
546 {
547
549 {
550 msp_timsData.get()
553 getMsRunId(), handler, ms_level);
554 }
555 }
556
557 catch(ExceptionInterrupted &)
558 {
559 qDebug() << "Reading of MS data interrupted by the user.";
560 }
561
562 // Now let the loading handler know that the loading of the data has ended.
563 // The handler might need this "signal" to perform additional tasks or to
564 // cleanup cruft.
565
566 // qDebug() << "Loading ended";
567 handler.loadingEnded();
568}
569
570
571std::size_t
576
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 setPrecursorNativeId(const QString &native_id)
Set the scan native id of the precursor ion.
void setDtInMilliSeconds(pappso_double rt)
Set the drift time in milliseconds.
void appendPrecursorIonData(const PrecursorIonData &precursor_ion_data)
const MassSpectrumId & getMassSpectrumId() const
Get the MassSpectrumId.
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.
pappso_double getRtInSeconds() const
Get 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
TimsDiaSlices * getTimsDiaSlicesPtr() const
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
bool isDiaRun() const
tells if this MS run is a DIA run
bool isDdaRun() const
tells if this MS run is a DDA run
TimsDdaPrecursors * getTimsDdaPrecursorsPtr() const
std::vector< TimsDdaPrecursors::SpectrumDescr > getSpectrumDescrListByFrameId(std::size_t frame_id) const
get a list of TimsDdaPrecursors::SpectrumDescr for a frame
void rawReaderSpectrumCollectionByMsLevel(const MsRunIdCstSPtr &msrun_id, SpectrumCollectionHandlerInterface &handler, unsigned int ms_level)
function to visit an MsRunReader and get each raw Spectrum in a spectrum collection handler by Ms Lev...
std::size_t getGlobalSliceIndexBeginByFrameId(std::size_t frame_id) const
std::size_t getLastMs1FrameIdByMs2FrameId(std::size_t frame_id) const
const std::map< std::size_t, MsMsWindowGroup * > & getMapFrame2WindowGroupPtr() const
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 getMobilityScan(std::size_t scanNum, std::size_t mz_index_merge_window, double mz_range_begin, double mz_range_end, quint32 &mz_minimum_index_out, quint32 &mz_maximum_index_out) const override
get a single mobility scan m/z + intensities
virtual void initialize() override
TimsMsRunReader(MsRunIdCstSPtr &msrun_id_csp)
virtual void readSpectrumCollection2(const MsRunReadConfig &config, SpectrumCollectionHandlerInterface &handler) override
virtual void readSpectrumCollection(SpectrumCollectionHandlerInterface &handler) override
function to visit an MsRunReader and get each Spectrum in a spectrum collection handler
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
virtual MassSpectrumCstSPtr massSpectrumCstSPtr(std::size_t spectrum_index) override
virtual MassSpectrumSPtr massSpectrumSPtr(std::size_t spectrum_index) override
get a MassSpectrumSPtr class given its spectrum index
virtual QualifiedMassSpectrum qualifiedMassSpectrum(std::size_t spectrum_index, bool want_binary_data=true) const override
get a QualifiedMassSpectrum class given its scan number
virtual std::size_t spectrumListSize() const override
get the totat number of spectrum conained in the MSrun data file
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
@ IsolationMzWidth
m/z isolation window width (left + right)
@ TimsFrameMzIndexBegin
Bruker's timsTOF mz index frame start range.
@ CollisionEnergy
Bruker's timsTOF collision energy.
@ TimsFrameMzIndexEnd
Bruker's timsTOF mz index frame end range.
@ BrukerPrecursorIndex
Bruker's timsTOF precursor index.
std::shared_ptr< MassSpectrum > MassSpectrumSPtr
std::shared_ptr< const TimsFrame > TimsFrameCstSPtr
Definition timsframe.h:43
std::size_t m_globalScanIndex
Definition timsdata.h:63
std::size_t m_scanCount
Definition timsdata.h:60
handle specific data for DDA MS runs
handle specific data for DIA MS runs
MSrun file reader for native Bruker TimsTOF raw data.