timsframebase.cpp

Go to the documentation of this file.

/**
 * \file pappsomspp/vendors/tims/timsframebase.cpp
 * \date 16/12/2019
 * \author Olivier Langella
 * \brief handle a single Bruker's TimsTof frame without binary 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 "timsframebase.h"
#include "pappsomspp/core/pappsoexception.h"
#include "pappsomspp/core/exception/exceptionoutofrange.h"
#include "mzcalibration/mzcalibrationmodel1.h"
#include <QDebug>
#include <QObject>
#include <cmath>
#include <algorithm>
 
namespace pappso
{
 
TimsFrameBase::TimsFrameBase(std::size_t frame_id, quint32 scanCount)
{
  qDebug() << frame_id;
  m_frameId = frame_id;
 
  m_scanCount = scanCount;
}
 
TimsFrameBase::TimsFrameBase([[maybe_unused]] const TimsFrameBase &other)
{
}
 
TimsFrameBase::~TimsFrameBase()
{
}
 
void
TimsFrameBase::setAcqDurationInMilliseconds(double acquisition_duration_ms)
{
  m_acqDurationInMilliseconds = acquisition_duration_ms;
}
 
 
void
TimsFrameBase::setMzCalibration(double T1_frame,
                                double T2_frame,
                                double digitizerTimebase,
                                double digitizerDelay,
                                double C0,
                                double C1,
                                double C2,
                                double C3,
                                double C4,
                                double T1_ref,
                                double T2_ref,
                                double dC1,
                                double dC2)
{
 
  /*  MzCalibrationModel1 mzCalibration(temperature_correction,
                                 digitizerTimebase,
                                 digitizerDelay,
                                 C0,
                                 C1,
                                 C2,
                                 C3,
                                 C4);
                                 */
  msp_mzCalibration = std::make_shared<MzCalibrationModel1>(T1_frame,
                                                            T2_frame,
                                                            digitizerTimebase,
                                                            digitizerDelay,
                                                            C0,
                                                            C1,
                                                            C2,
                                                            C3,
                                                            C4,
                                                            T1_ref,
                                                            T2_ref,
                                                            dC1,
                                                            dC2);
}
 
bool
TimsFrameBase::checkScanNum(std::size_t scanNum) const
{
  if(scanNum >= m_scanCount)
    {
      throw pappso::ExceptionOutOfRange(
        QObject::tr("Invalid scan number : scanNum %1  > m_scanNumber %2")
          .arg(scanNum)
          .arg(m_scanCount));
    }
 
  return true;
}
 
std::size_t
TimsFrameBase::getScanPeakCount(std::size_t scanNum) const
{
  throw PappsoException(
    QObject::tr("ERROR unable to get number of peaks in TimsFrameBase for scan number %1")
      .arg(scanNum));
}
 
std::size_t
TimsFrameBase::getTotalNumberOfScans() const
{
  return m_scanCount;
}
 
MassSpectrumSPtr
TimsFrameBase::getMassSpectrumSPtr(std::size_t scanNum) const
{
  throw PappsoException(
    QObject::tr("ERROR unable to getMassSpectrumSPtr in TimsFrameBase for scan number %1")
      .arg(scanNum));
}
 
 
pappso::MassSpectrumCstSPtr
TimsFrameBase::getMassSpectrumCstSPtr(std::size_t scanNum) const
{
  // qDebug();
 
  return getMassSpectrumSPtr(scanNum);
}
Trace
TimsFrameBase::cumulateScansToTrace(std::size_t scanNumBegin, std::size_t scanNumEnd) const
{
  throw PappsoException(QObject::tr("ERROR unable to cumulateScanToTrace in TimsFrameBase for scan "
                                    "number begin %1 end %2")
                          .arg(scanNumBegin)
                          .arg(scanNumEnd));
}
 
Trace
TimsFrameBase::combineScansToTraceWithDowngradedMzResolution(
  std::size_t mzindex_merge_window [[maybe_unused]],
  std::size_t scanNumBegin [[maybe_unused]],
  std::size_t scanNumEnd [[maybe_unused]],
  quint32 &minimum_tof_index_out [[maybe_unused]],
  quint32 &maximum_tof_index_out [[maybe_unused]]) const
{
  throw PappsoException(
    QObject::tr("Non implemented function %1 %2 %3").arg(__FILE__).arg(__FUNCTION__).arg(__LINE__));
}
 
Trace
TimsFrameBase::combineScansToTraceWithDowngradedMzResolution2(
  std::size_t mz_index_merge_window [[maybe_unused]],
  double mz_range_begin [[maybe_unused]],
  double mz_range_end [[maybe_unused]],
  std::size_t mobility_scan_begin [[maybe_unused]],
  std::size_t mobility_scan_end [[maybe_unused]],
  quint32 &mz_minimum_index_out [[maybe_unused]],
  quint32 &mz_maximum_index_out [[maybe_unused]]) const
{
  throw PappsoException(
    QObject::tr("Non implemented function %1 %2 %3").arg(__FILE__).arg(__FUNCTION__).arg(__LINE__));
}
 
Trace
TimsFrameBase::getMobilityScan(std::size_t scanNum [[maybe_unused]],
                               std::size_t mz_index_merge_window [[maybe_unused]],
                               double mz_range_begin [[maybe_unused]],
                               double mz_range_end [[maybe_unused]],
                               quint32 &mz_minimum_index_out [[maybe_unused]],
                               quint32 &mz_maximum_index_out [[maybe_unused]]) const
{
  throw PappsoException(
    QObject::tr("Non implemented function %1 %2 %3").arg(__FILE__).arg(__FUNCTION__).arg(__LINE__));
}
 
void
TimsFrameBase::combineScansInTofIndexIntensityMap(TimsDataFastMap &rawSpectrum [[maybe_unused]],
                                                  std::size_t scanNumBegin,
                                                  std::size_t scanNumEnd) const
{
  throw PappsoException(
    QObject::tr("ERROR unable to cumulateScansInRawMap in TimsFrameBase for scan "
                "number begin %1 end %2")
      .arg(scanNumBegin)
      .arg(scanNumEnd));
}
 
 
void
TimsFrameBase::combineScansInTofIndexIntensityMap(TimsDataFastMap &rawSpectrum [[maybe_unused]],
                                                  std::size_t scanNumBegin,
                                                  std::size_t scanNumEnd,
                                                  quint32 tof_index_begin,
                                                  quint32 tof_index_end) const
{
  throw PappsoException(
    QObject::tr("ERROR unable to cumulateScansInRawMap in TimsFrameBase for scan "
                "number begin %1 end %2, tof index begin %3, tof index end %4")
      .arg(scanNumBegin)
      .arg(scanNumEnd)
      .arg(tof_index_begin)
      .arg(tof_index_end));
}
 
quint64
TimsFrameBase::cumulateScanIntensities(std::size_t scanNum) const
{
  throw PappsoException(
    QObject::tr("ERROR unable to cumulateSingleScanIntensities in TimsFrameBase for scan "
                "number %1.")
      .arg(scanNum));
 
  return 0;
}
 
 
quint64
TimsFrameBase::cumulateScanRangeIntensities(std::size_t scanNumBegin, std::size_t scanNumEnd) const
{
  throw PappsoException(
    QObject::tr("ERROR unable to cumulateScansInRawMap in TimsFrameBase for scan "
                "number begin %1 end %2")
      .arg(scanNumBegin)
      .arg(scanNumEnd));
 
  return 0;
}
 
void
TimsFrameBase::setRtInSeconds(double time)
{
  m_rtInSeconds = time;
}
 
void
TimsFrameBase::setMsMsType(quint8 type)
{
 
  qDebug() << " m_msMsType=" << type;
  m_msMsType = type;
}
 
unsigned int
TimsFrameBase::getMsLevel() const
{
  if(m_msMsType == 0)
    return 1;
  return 2;
}
 
double
TimsFrameBase::getRtInSeconds() const
{
  return m_rtInSeconds;
}
 
std::size_t
TimsFrameBase::getId() const
{
  return m_frameId;
}
void
TimsFrameBase::setTimsCalibration(int tims_model_type,
                                  double C0,
                                  double C1,
                                  double C2,
                                  double C3,
                                  double C4,
                                  [[maybe_unused]] double C5,
                                  double C6,
                                  double C7,
                                  double C8,
                                  double C9)
{
  if(tims_model_type != 2)
    {
      throw pappso::PappsoException(
        QObject::tr("ERROR in TimsFrame::setTimsCalibration tims_model_type != 2"));
    }
  m_timsDvStart = C2; // C2 from TimsCalibration
  m_timsTtrans  = C4; // C4 from TimsCalibration
  m_timsNdelay  = C0; // C0 from TimsCalibration
  m_timsVmin    = C8; // C8 from TimsCalibration
  m_timsVmax    = C9; // C9 from TimsCalibration
  m_timsC6      = C6;
  m_timsC7      = C7;
 
 
  m_timsSlope = (C3 - m_timsDvStart) / C1; //  //C3 from TimsCalibration // C2 from
                                           //  TimsCalibration // C1 from TimsCalibration
}
double
TimsFrameBase::getVoltageTransformation(std::size_t scanNum) const
{
  double v = m_timsDvStart + m_timsSlope * ((double)scanNum - m_timsTtrans - m_timsNdelay);
 
  if(v < m_timsVmin)
    {
      throw pappso::PappsoException(
        QObject::tr("ERROR in TimsFrame::getVoltageTransformation invalid tims "
                    "calibration, v < m_timsVmin %1 < %2")
          .arg(v)
          .arg(m_timsVmin));
    }
 
 
  if(v > m_timsVmax)
    {
      throw pappso::PappsoException(
        QObject::tr("ERROR in TimsFrame::getVoltageTransformation invalid tims "
                    "calibration, v > m_timsVmax  %1 > %2")
          .arg(v)
          .arg(m_timsVmax));
    }
  return v;
}
double
TimsFrameBase::getDriftTimeInMilliseconds(std::size_t scanNum) const
{
  return (m_acqDurationInMilliseconds / (double)m_scanCount) * ((double)scanNum);
}
 
double
TimsFrameBase::getOneOverK0Transformation(std::size_t scanNum) const
{
  return 1 / (m_timsC6 + (m_timsC7 / getVoltageTransformation(scanNum)));
}
 
 
std::size_t
TimsFrameBase::getScanIndexFromOneOverK0(double one_over_k0) const
{
  double temp = 1 / one_over_k0;
  temp        = temp - m_timsC6;
  temp        = temp / m_timsC7;
  temp        = 1 / temp;
  temp        = temp - m_timsDvStart;
  temp        = temp / m_timsSlope + m_timsTtrans + m_timsNdelay;
  return (std::size_t)std::round(temp);
}
 
bool
TimsFrameBase::hasSameCalibrationData(const TimsFrameBase &other) const
{
  if((m_timsDvStart == other.m_timsDvStart) && (m_timsTtrans == other.m_timsTtrans) &&
     (m_timsNdelay == other.m_timsNdelay) && (m_timsVmin == other.m_timsVmin) &&
     (m_timsVmax == other.m_timsVmax) && (m_timsC6 == other.m_timsC6) &&
     (m_timsC7 == other.m_timsC7) && (m_timsSlope == other.m_timsSlope))
    {
      return true;
    }
  return false;
}
 
 
pappso::Trace
TimsFrameBase::getTraceFromTofIndexIntensityMap(TimsDataFastMap &accumulated_scans) const
{
  qDebug();
  // qDebug();
  // add flanking peaks
  pappso::Trace local_trace;
 
  MzCalibrationInterface *mz_calibration_p = getMzCalibrationInterfaceSPtr().get();
 
 
  DataPoint element;
  for(quint32 tof_index : accumulated_scans.getTofIndexList())
    {
      // intensity normalization
      element.y =
        ((double)accumulated_scans.readIntensity(tof_index)) * 100.0 / m_acqDurationInMilliseconds;
 
      // mz calibration
      element.x = mz_calibration_p->getMzFromTofIndex(tof_index);
 
      local_trace.push_back(element);
    }
  local_trace.sortX();
 
  qDebug();
  // qDebug();
  return local_trace;
}
 
const MzCalibrationInterfaceSPtr &
TimsFrameBase::getMzCalibrationInterfaceSPtr() const
{
  if(msp_mzCalibration == nullptr)
    {
 
      throw pappso::PappsoException(QObject::tr("ERROR in %1, %2, %3 msp_mzCalibration is null")
                                      .arg(__FILE__)
                                      .arg(__FUNCTION__)
                                      .arg(__LINE__));
    }
  return msp_mzCalibration;
}
 
void
TimsFrameBase::setMzCalibrationInterfaceSPtr(MzCalibrationInterfaceSPtr mzCalibration)
{
 
  if(mzCalibration == nullptr)
    {
 
      throw pappso::PappsoException(QObject::tr("ERROR in %1, %2, %3 msp_mzCalibration is null")
                                      .arg(__FILE__)
                                      .arg(__FUNCTION__)
                                      .arg(__LINE__));
    }
  msp_mzCalibration = mzCalibration;
}
 
 
quint32
TimsFrameBase::getMaximumRawMassIndex() const
{
  quint32 max_value = 0;
  for(quint32 i = 0; i < m_scanCount; i++)
    {
      qDebug() << "m_scanNumber=" << m_scanCount << " i=" << i;
      std::vector<quint32> index_list = getScanTofIndexList(i);
      auto it                         = std::max_element(index_list.begin(), index_list.end());
      if(it != index_list.end())
        {
          max_value = std::max(max_value, *it);
        }
    }
  return max_value;
}
 
std::vector<quint32>
TimsFrameBase::getScanTofIndexList(std::size_t scanNum) const
{
  throw PappsoException(
    QObject::tr("ERROR unable to getScanIndexList in TimsFrameBase for scan number %1")
      .arg(scanNum));
}
 
 
std::vector<quint32>
TimsFrameBase::getScanIntensityList(std::size_t scanNum) const
{
  throw PappsoException(
    QObject::tr("ERROR unable to getScanIntensities in TimsFrameBase for scan number %1")
      .arg(scanNum));
}
 
Trace
TimsFrameBase::getIonMobilityTraceByTofIndexRange(std::size_t tof_index_begin,
                                                  std::size_t tof_index_end,
                                                  Enums::XicExtractMethod method,
                                                  std::size_t scan_index_begin,
                                                  std::size_t scan_index_end) const
{
  Trace im_trace;
  DataPoint data_point;
  if(scan_index_end > m_scanCount)
    scan_index_end = m_scanCount;
  for(quint32 i = scan_index_begin; i < (scan_index_end + 1); i++)
    {
      data_point.x = i;
      data_point.y = 0;
      qDebug() << "m_scanNumber=" << m_scanCount << " i=" << i;
      std::vector<quint32> index_list = getScanTofIndexList(i);
      auto it_lower =
        std::find_if(index_list.begin(), index_list.end(), [tof_index_begin](quint32 to_compare) {
          if(to_compare < tof_index_begin)
            {
              return false;
            }
          return true;
        });
 
 
      if(it_lower == index_list.end())
        {
        }
      else
        {
 
 
          auto it_upper =
            std::find_if(index_list.begin(), index_list.end(), [tof_index_end](quint32 to_compare) {
              if(tof_index_end >= to_compare)
                {
                  return false;
                }
              return true;
            });
          std::vector<quint32> intensity_list = getScanIntensityList(i);
          for(int j = std::distance(index_list.begin(), it_lower);
              j < std::distance(index_list.begin(), it_upper);
              j++)
            {
              if(method == Enums::XicExtractMethod::sum)
                {
                  data_point.y += intensity_list[j];
                }
              else
                {
                  data_point.y = std::max((double)intensity_list[j], data_point.y);
                }
            }
        }
      im_trace.push_back(data_point);
    }
  qDebug();
  return im_trace;
}
 
 
std::vector<TimsFrameBase::TofIndexIntensityPair> &
TimsFrameBase::downgradeResolutionOfTofIndexIntensityPairList(
  std::size_t mzindex_merge_window,
  std::vector<TimsFrameBase::TofIndexIntensityPair> &rawSpectrum) const
{
 
  qDebug() << rawSpectrum.size();
  std::vector<TimsFrameBase::TofIndexIntensityPair> new_spectrum;
 
  TimsFrameBase::TofIndexIntensityPair current_point;
  current_point.intensity_index = 0;
  current_point.tof_index       = 0;
  for(auto &pair_mz_intensity : rawSpectrum)
    {
      quint32 mzkey = (pair_mz_intensity.tof_index / mzindex_merge_window);
      mzkey         = (mzkey * mzindex_merge_window) + (mzindex_merge_window / 2);
 
      if(current_point.tof_index != mzkey)
        {
          if(current_point.tof_index > 0)
            {
              new_spectrum.push_back(current_point);
            }
 
          current_point.intensity_index = pair_mz_intensity.intensity_index;
          current_point.tof_index       = mzkey;
        }
      else
        {
          current_point.intensity_index += pair_mz_intensity.intensity_index;
        }
    }
 
  if(current_point.tof_index > 0)
    {
      new_spectrum.push_back(current_point);
    }
  rawSpectrum = new_spectrum;
  qDebug() << rawSpectrum.size();
  return rawSpectrum;
}
// namespace pappso
} // namespace pappso