timsframesmsrunreader.cpp

Go to the documentation of this file.

/**
 * \file pappsomspp/msrun/private/timsmsrunreader.h
 * \date 05/09/2019
 * \author Olivier Langella
 * \brief MSrun file reader for native Bruker TimsTOF raw data
 */
 
/*******************************************************************************
 * Copyright (c) 2019 Olivier Langella <Olivier.Langella@u-psud.fr>.
 *
 * This file is part of the PAPPSOms++ library.
 *
 *     PAPPSOms++ is free software: you can redistribute it and/or modify
 *     it under the terms of the GNU General Public License as published by
 *     the Free Software Foundation, either version 3 of the License, or
 *     (at your option) any later version.
 *
 *     PAPPSOms++ is distributed in the hope that it will be useful,
 *     but WITHOUT ANY WARRANTY; without even the implied warranty of
 *     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *     GNU General Public License for more details.
 *
 *     You should have received a copy of the GNU General Public License
 *     along with PAPPSOms++.  If not, see <http://www.gnu.org/licenses/>.
 *
 ******************************************************************************/
 
#include "timsframesmsrunreader.h"
#include "pappsomspp/core/exception/exceptionnotimplemented.h"
#include "pappsomspp/core/exception/exceptioninterrupted.h"
#include "pappsomspp/core/exception/exceptionnotpossible.h"
#include "../../vendors/tims/timsddaprecursors.h"
#include <QDebug>
 
using namespace pappso;
 
TimsFramesMsRunReader::TimsFramesMsRunReader(MsRunIdCstSPtr &msrun_id_csp)
  : TimsMsRunReaderBase(msrun_id_csp)
{
}
 
TimsFramesMsRunReader::~TimsFramesMsRunReader()
{
  msp_timsData = nullptr;
}
 
 
pappso::MassSpectrumSPtr
TimsFramesMsRunReader::massSpectrumSPtr([[maybe_unused]] std::size_t spectrum_index)
{
  throw ExceptionNotImplemented(
    QObject::tr("Not yet implemented in TimsFramesMsRunReader %1.\n").arg(__LINE__));
 
  return pappso::MassSpectrumSPtr();
}
 
 
pappso::MassSpectrumCstSPtr
TimsFramesMsRunReader::massSpectrumCstSPtr(std::size_t spectrum_index)
{
  return msp_timsData->getMassSpectrumCstSPtrByGlobalScanIndex(spectrum_index);
}
 
 
QualifiedMassSpectrum
TimsFramesMsRunReader::qualifiedMassSpectrum(std::size_t spectrum_index,
                                             bool want_binary_data) const
{
 
  QualifiedMassSpectrum mass_spectrum;
 
  msp_timsData->getQualifiedMassSpectrumByGlobalScanIndex(
    getMsRunId(), mass_spectrum, spectrum_index, want_binary_data);
  return mass_spectrum;
}
 
 
void
TimsFramesMsRunReader::readSpectrumCollection(SpectrumCollectionHandlerInterface &handler)
{
  // qDebug() << "Reading the spectrum collection with no specific
  // configuration.";
  MsRunReadConfig config;
  readSpectrumCollection2(config, handler);
}
 
 
void
TimsFramesMsRunReader::readSpectrumCollection2(const MsRunReadConfig &config,
                                               SpectrumCollectionHandlerInterface &handler)
{
  // qDebug().noquote() << "Reading the spectrum collection with this "
  //                       "specific configuration:"
  //                    << config.toString();
 
  std::vector<std::size_t> subset_of_tims_frame_ids;
 
 
  bool asked_ion_mobility_scan_num_range = false;
 
  quint32 mobility_scan_num_range_begin = std::numeric_limits<quint32>::max();
  quint32 mobility_scan_num_range_end   = std::numeric_limits<quint32>::max();
  quint32 mobility_scan_num_range_width = std::numeric_limits<quint32>::max();
 
  double mobility_one_over_k0             = std::numeric_limits<double>::max();
  double mobility_one_over_k0_range_begin = std::numeric_limits<double>::max();
  double mobility_one_over_k0_range_end   = std::numeric_limits<double>::max();
 
  if(!config.getParameterValue(MsRunReadConfigParameter::TimsFrameIonMobScanIndexBegin).isNull() &&
     !config.getParameterValue(MsRunReadConfigParameter::TimsFrameIonMobScanIndexEnd).isNull())
    {
      mobility_scan_num_range_begin =
        config.getParameterValue(MsRunReadConfigParameter::TimsFrameIonMobScanIndexBegin).toUInt();
      mobility_scan_num_range_end =
        config.getParameterValue(MsRunReadConfigParameter::TimsFrameIonMobScanIndexEnd).toUInt();
 
      // We need the range width below.
      mobility_scan_num_range_width =
        mobility_scan_num_range_end + 1 - mobility_scan_num_range_begin;
 
      asked_ion_mobility_scan_num_range = true;
 
      // Be sure to check in the frames loop below that the user might
      // have asked for an ion mobility range but on the basis of the 1/K0 unit.
    }
 
  const std::vector<FrameIdDescr> &frame_id_descr_list = msp_timsData->getFrameIdDescrList();
 
  // Just for the feedback to the user.
  std::size_t scan_count = 0;
 
  for(auto const &frame_record : msp_timsData->getTimsFrameRecordList())
    {
      if(handler.shouldStop())
        {
          // qDebug() << "The operation was cancelled. Breaking the loop.";
          throw ExceptionInterrupted(QObject::tr("Reading timsTOF data cancelled by the user."));
        }
 
      if(frame_record.frame_id == 0)
        continue;
 
      if(!config.acceptRetentionTimeInSeconds(frame_record.frame_time))
        continue;
 
      std::size_t ms_level = 2;
      if(frame_record.msms_type == 0)
        ms_level = 1;
 
      if(!config.acceptMsLevel(ms_level))
        continue;
 
      // FIXME: this might be the place where to actually get the list of
      // precursors' m/z and charge values.
 
      subset_of_tims_frame_ids.push_back(frame_record.frame_id);
 
      if(mobility_scan_num_range_width != std::numeric_limits<int>::max())
        {
          scan_count += mobility_scan_num_range_width;
        }
      else
        {
          scan_count += frame_id_descr_list[frame_record.frame_id].m_scanCount;
        }
    }
 
  // At this point, we have a subset of frame records.
  std::size_t frame_count = subset_of_tims_frame_ids.size();
  qDebug() << "The number of retained RT range- and MS level-matching frames : " << frame_count;
 
  // Inform the handler of the spectrum list so that it can handle feedback to
  // the user.
 
  // FIXME:
  // Either we document the number of frames (because we assume we will
  // flatten them all)...
  handler.spectrumListHasSize(frame_count);
  // Or we document the number of actual scans because we might not flatten
  // all the frames.
  // handler.spectrumListHasSize(scan_count);
 
  // Check for m/z range selection
  bool asked_mz_range   = false;
  double mz_range_begin = -1;
  double mz_range_end   = -1;
 
  if(!config.getParameterValue(MsRunReadConfigParameter::MzRangeBegin).isNull())
    {
      asked_mz_range = true;
 
      mz_range_begin = config.getParameterValue(MsRunReadConfigParameter::MzRangeBegin).toDouble();
      mz_range_end   = config.getParameterValue(MsRunReadConfigParameter::MzRangeEnd).toDouble();
 
      // qDebug() << "The m/z range asked is: " << mz_range_begin << "--"
      //          << mz_range_end;
    }
 
  // Check for m/z resolution downgrading
  // The idea is that we merge a number of mz indices into a single index,
  // which is essentially an increase of the m/z bin size, and therefore
  // of the resolution/definition of the mass spectrum.
  std::size_t mz_index_merge_window = 0;
  if(!config.getParameterValue(MsRunReadConfigParameter::TimsFrameMzIndexMergeWindow).isNull())
    {
      mz_index_merge_window =
        config.getParameterValue(MsRunReadConfigParameter::TimsFrameMzIndexMergeWindow).toUInt();
 
      // qDebug() << "mz_index_merge_window=" << mz_index_merge_window;
    }
 
  // The scan index is the index of the scan in the *whole* mass data file, it
  // is a sequential number of scans over all the frames.
  std::size_t scan_index = 0; // iterate in each spectrum
 
  for(std::size_t tims_frame_id : subset_of_tims_frame_ids)
    {
      if(handler.shouldStop())
        {
          // qDebug() << "The operation was cancelled. Breaking the loop.";
          throw ExceptionInterrupted(QObject::tr("Reading timsTOF data cancelled by the user."));
        }
 
      // qDebug() << "tims_frame_id=" << tims_frame_id;
 
      const FrameIdDescr &current_frame_record = frame_id_descr_list[tims_frame_id];
      // qDebug() << "tims_frame_id=" << tims_frame_id;
 
      // FIXME: from an outsider point of view, cumulated size does not
      // convey the notion of sequential scan number.
 
      scan_index = current_frame_record.m_globalScanIndex;
 
      TimsFrameCstSPtr tims_frame_csp = msp_timsData->getTimsFrameCstSPtrCached(tims_frame_id);
 
      // If the user wants to select 1/Ko values in a given range, we need to
      // compute the ion mobility scan value starting from that 1/Ko value in
      // *each* frame. Note that the computed mobility_scan_num_begin and
      // mobility_scan_num_end would override thoses possibly set with
      // TimsFramesMsRunReader_mobility_index_begin/end above.
 
      if(!config.getParameterValue(MsRunReadConfigParameter::TimsFrameIonMobOneOverK0Begin)
            .isNull() &&
         !config.getParameterValue(MsRunReadConfigParameter::TimsFrameIonMobOneOverK0End).isNull())
        {
          mobility_one_over_k0_range_begin =
            config.getParameterValue(MsRunReadConfigParameter::TimsFrameIonMobOneOverK0Begin)
              .toDouble();
 
          mobility_one_over_k0_range_end =
            config.getParameterValue(MsRunReadConfigParameter::TimsFrameIonMobOneOverK0End)
              .toDouble();
 
          mobility_scan_num_range_begin =
            tims_frame_csp.get()->getScanIndexFromOneOverK0(mobility_one_over_k0_range_begin);
 
          mobility_scan_num_range_end =
            tims_frame_csp.get()->getScanIndexFromOneOverK0(mobility_one_over_k0_range_end);
 
          asked_ion_mobility_scan_num_range = true;
        }
 
      // Now that we know if the user has asked for an ion mobility range,
      // either using scan indices or 1/K0 values, we need to double check the
      // range borders.
      quint32 count_of_mobility_scans = tims_frame_csp->getTotalNumberOfScans();
 
      if(asked_ion_mobility_scan_num_range)
        {
          if(mobility_scan_num_range_end > (count_of_mobility_scans - 1))
            {
              mobility_scan_num_range_end = count_of_mobility_scans - 1;
            }
        }
      else
        {
          mobility_scan_num_range_begin = 0;
          mobility_scan_num_range_end   = count_of_mobility_scans - 1;
        }
 
      // Now that we know the mobility index range, if we did not set the
      // mobility one over K0 because that was not the unit used by
      // the caller, then we can compute these values and set them
      // later to the qualified mass spectrum parameters.
      if(mobility_one_over_k0_range_begin == std::numeric_limits<double>::max())
        mobility_one_over_k0_range_begin =
          tims_frame_csp->getOneOverK0Transformation(mobility_scan_num_range_begin);
      if(mobility_one_over_k0_range_end == std::numeric_limits<double>::max())
        mobility_one_over_k0_range_end =
          tims_frame_csp->getOneOverK0Transformation(mobility_scan_num_range_end);
 
      mobility_scan_num_range_width =
        mobility_scan_num_range_end + 1 - mobility_scan_num_range_begin;
 
      // We want to provide the inverse mobility for the scan that sits in the
      // middle of the defined range or the whole range if none is defined..
      mobility_one_over_k0 = tims_frame_csp.get()->getScanIndexFromOneOverK0(
        mobility_scan_num_range_begin + (mobility_scan_num_range_width / 2));
 
      // Now, with or without the peak list, we have to craft a qualified mass
      // spectrum that will hold all the data about the data in it.
      QualifiedMassSpectrum qualified_mass_spectrum;
 
      MassSpectrumId spectrum_id;
 
      spectrum_id.setSpectrumIndex(tims_frame_id);
      spectrum_id.setMsRunId(getMsRunId());
 
      // Can be modified to add bits that might help our case
      spectrum_id.setNativeId(QString("frame_id=%1 global_scan_index=%2 im_scan_range_begin=%3 "
                                      "im_scan_range_end=%4")
                                .arg(tims_frame_id)
                                .arg(scan_index)
                                .arg(mobility_scan_num_range_begin)
                                .arg(mobility_scan_num_range_end));
 
      qualified_mass_spectrum.setMassSpectrumId(spectrum_id);
 
      // We want to document the retention time!
 
      qualified_mass_spectrum.setRtInSeconds(tims_frame_csp.get()->getRtInSeconds());
 
      // We do want to document the ms level of the spectrum and possibly
      // the precursor's m/z and charge.
      unsigned int frame_ms_level = tims_frame_csp.get()->getMsLevel();
      qualified_mass_spectrum.setMsLevel(frame_ms_level);
 
 
      // Arrival time at half the range.
 
      qualified_mass_spectrum.setDtInMilliSeconds(tims_frame_csp.get()->getDriftTimeInMilliseconds(
        mobility_scan_num_range_begin + (mobility_scan_num_range_width / 2)));
 
      // 1/K0
      qDebug() << "mobility_one_over_k0:" << mobility_one_over_k0
               << "mobility_one_over_k0_range_begin:" << mobility_one_over_k0_range_begin
               << "mobility_one_over_k0_range_end" << mobility_one_over_k0_range_end;
 
      if(mobility_one_over_k0 == std::numeric_limits<double>::max() ||
         mobility_one_over_k0_range_begin == std::numeric_limits<double>::max() ||
         mobility_one_over_k0_range_end == std::numeric_limits<double>::max())
        throw(
          ExceptionNotPossible("Not possible that mobility_one_over_k0 and its "
                               "range are undefined."));
 
      qualified_mass_spectrum.setParameterValue(QualifiedMassSpectrumParameter::IonMobOneOverK0,
                                                mobility_one_over_k0);
      qualified_mass_spectrum.setParameterValue(
        QualifiedMassSpectrumParameter::IonMobOneOverK0Begin, mobility_one_over_k0_range_begin);
      qualified_mass_spectrum.setParameterValue(QualifiedMassSpectrumParameter::IonMobOneOverK0End,
                                                mobility_one_over_k0_range_end);
 
      // qDebug() << "mobility_scan_num_range_begin:"
      //          << mobility_scan_num_range_begin
      //          << "mobility_scan_num_range_end:" <<
      //          mobility_scan_num_range_end;
 
      if(mobility_scan_num_range_begin == std::numeric_limits<quint32>::max() ||
         mobility_scan_num_range_end == std::numeric_limits<quint32>::max())
        throw(
          ExceptionNotPossible("Not possible that mobility_scan_num_range values are undefined."));
 
      qualified_mass_spectrum.setParameterValue(QualifiedMassSpectrumParameter::TimsIonMobScanIndex,
                                                mobility_scan_num_range_begin +
                                                  (mobility_scan_num_range_width / 2));
      qualified_mass_spectrum.setParameterValue(
        QualifiedMassSpectrumParameter::TimsFrameIonMobScanIndexBegin,
        mobility_scan_num_range_begin);
      qualified_mass_spectrum.setParameterValue(
        QualifiedMassSpectrumParameter::TimsFrameIonMobScanIndexEnd, mobility_scan_num_range_end);
 
      qualified_mass_spectrum.setParameterValue(
        QualifiedMassSpectrumParameter::TimsFrameScansCount,
        static_cast<qlonglong>(tims_frame_csp->getTotalNumberOfScans()));
 
 
      Trace trace;
 
      if(config.needPeakList())
        {
          // Provide these two variables for the function below to fill in the
          // values. We will need them later.
          quint32 min_mz_index_out = 0;
          quint32 max_mz_index_out = 0;
 
          if(asked_mz_range)
            {
              trace = tims_frame_csp->combineScansToTraceWithDowngradedMzResolution2(
                mz_index_merge_window,
                mz_range_begin,
                mz_range_end,
                mobility_scan_num_range_begin,
                mobility_scan_num_range_end,
                min_mz_index_out,
                max_mz_index_out);
            }
          else
            {
              trace = tims_frame_csp->combineScansToTraceWithDowngradedMzResolution(
                mz_index_merge_window,
                mobility_scan_num_range_begin,
                mobility_scan_num_range_end,
                min_mz_index_out,
                max_mz_index_out);
            }
 
          // qDebug() << "Got min_mz_index_out:" << min_mz_index_out;
          // qDebug() << "Got max_mz_index_out:" << max_mz_index_out;
 
          qualified_mass_spectrum.setParameterValue(
            QualifiedMassSpectrumParameter::TimsFrameMzIndexBegin, min_mz_index_out);
          qualified_mass_spectrum.setParameterValue(
            QualifiedMassSpectrumParameter::TimsFrameMzIndexEnd, max_mz_index_out);
 
          qualified_mass_spectrum.setMassSpectrumSPtr(std::make_shared<MassSpectrum>(trace));
          qualified_mass_spectrum.setEmptyMassSpectrum(false);
        }
      else
        {
          qualified_mass_spectrum.setEmptyMassSpectrum(true);
        }
 
      handler.setQualifiedMassSpectrum(qualified_mass_spectrum);
    }
}
 
 
void
TimsFramesMsRunReader::readSpectrumCollectionByMsLevel(
  [[maybe_unused]] SpectrumCollectionHandlerInterface &handler,
  [[maybe_unused]] unsigned int ms_level)
{
  qDebug();
}
 
 
std::size_t
TimsFramesMsRunReader::spectrumListSize() const
{
  return msp_timsData->getTotalScanCount();
}
 
 
Trace
TimsFramesMsRunReader::computeTicChromatogram()
{
 
  // We want to compute the TIC chromatogram, not load the chromatogram that
  // is located in the SQL database.
  //
  // For this, we need to iterated into the frames and ask for MS1 spectra
  // only. msp_timsData has that information:
  //
  // std::vector<FrameIdDescr> m_frameIdDescrList;
  //
  // and
 
  // struct FrameIdDescr
  // {
  //   std::size_t m_frameId;   // frame id
  //   std::size_t m_size;      // frame size (number of TOF scans in frame)
  //   std::size_t m_cumulSize; // cumulative size
  // };
 
  Trace tic_chromatogram;
 
  const std::vector<FrameIdDescr> frame_descr_list = msp_timsData->getFrameIdDescrList();
 
  for(FrameIdDescr frame_id_descr : frame_descr_list)
    {
      TimsFrameCstSPtr tims_frame_csp =
        msp_timsData->getTimsFrameCstSPtrCached(frame_id_descr.m_frameId);
      std::size_t scan_begin = 0;
      std::size_t scan_end   = tims_frame_csp->getTotalNumberOfScans() - 1;
 
      // By convention, a TIC chromatogram is only performed using MS1
      // spectra.
      if(tims_frame_csp->getMsLevel() == 1)
        {
 
          // Retention times are in seconds in the Bruker world.
          double rt = tims_frame_csp->getRtInSeconds();
 
          tic_chromatogram.append(
            DataPoint(rt, tims_frame_csp->cumulateScanRangeIntensities(scan_begin, scan_end)));
        }
      else
        continue;
    }
 
  return tic_chromatogram;
}
 
 
std::size_t
TimsFramesMsRunReader::spectrumStringIdentifier2SpectrumIndex(const QString &spectrum_identifier
                                                              [[maybe_unused]])
{
  throw pappso::ExceptionNotImplemented(
    QObject::tr("%1 %2 %3 not implemented").arg(__FILE__).arg(__FUNCTION__).arg(__LINE__));
}