timsdiaslices.cpp

Go to the documentation of this file.

/**
 * \file pappsomspp/vendors/tims/timsdiaslices.cpp
 * \date 02/07/2024
 * \brief handle specific data for DIA MS runs
 */
 
/*******************************************************************************
 * Copyright (c) 2024 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 "timsdiaslices.h"
#include "pappsomspp/core/exception/exceptionnotfound.h"
 
#include <iterator>
#include <QSqlError>
 
namespace pappso
{
 
bool
TimsDiaSlices::MsMsWindow::isMzInRange(double mz) const
{
  double half_window = IsolationWidth / 2;
  if(mz < (IsolationMz - half_window))
    return false;
  if(mz > (IsolationMz + half_window))
    return false;
  return true;
}
 
TimsDiaSlices::MsMsWindowGroupList::~MsMsWindowGroupList()
{
  for(MsMsWindowGroup *p_group : *this)
    {
      delete p_group;
    }
}
 
void
TimsDiaSlices::MsMsWindowGroupList::addInGroup(const MsMsWindow &window)
{
  std::size_t group_id = window.WindowGroup;
  auto it_in_list      = std::find_if(begin(), end(), [group_id](const MsMsWindowGroup *group) {
    if(group->WindowGroup == group_id)
      return true;
    return false;
  });
  if(it_in_list == end())
    {
 
      qDebug();
      MsMsWindowGroup *new_group = new MsMsWindowGroup();
      new_group->WindowGroup     = group_id;
      new_group->push_back(window);
      push_back(new_group);
    }
  else
    {
      qDebug();
      (*it_in_list)->push_back(window);
 
      std::sort(
        (*it_in_list)->begin(), (*it_in_list)->end(), [](const MsMsWindow &a, const MsMsWindow &b) {
          return (a.ScanNumBegin < b.ScanNumBegin);
        });
      std::size_t i = 0;
      for(MsMsWindow &window : *(*it_in_list))
        {
          window.SliceIndex = i;
          i++;
        }
    }
  qDebug();
}
 
TimsDiaSlices::MsMsWindowGroup *
TimsDiaSlices::MsMsWindowGroupList::getWindowGroupPtrByGroupId(std::size_t window_group_id) const
{
 
  auto it_in_list = std::find_if(begin(), end(), [window_group_id](const MsMsWindowGroup *group) {
    if(group->WindowGroup == window_group_id)
      return true;
    return false;
  });
 
  if(it_in_list == end())
    {
      throw pappso::ExceptionNotFound(
        QObject::tr("MsMsWindowGroup %1 not found").arg(window_group_id));
    }
  return *it_in_list;
}
 
TimsDiaSlices::TimsDiaSlices(QSqlQuery &query, TimsData *tims_data_origin)
  : mp_timsDataOrigin(tims_data_origin)
{
 
  /*
CREATE TABLE IF NOT EXISTS "DiaFrameMsMsWindows" (
"WindowGroup"   INTEGER NOT NULL,
"ScanNumBegin"  INTEGER NOT NULL,
"ScanNumEnd"    INTEGER NOT NULL,
"IsolationMz"   REAL NOT NULL,
"IsolationWidth"    REAL NOT NULL,
"CollisionEnergy"   REAL NOT NULL,
FOREIGN KEY("WindowGroup") REFERENCES "FrameMsMsWindowGroups"("Id"),
PRIMARY KEY("WindowGroup","ScanNumBegin")
) WITHOUT ROWID;
*/
 
  if(!query.exec("SELECT WindowGroup, ScanNumBegin, ScanNumEnd, IsolationMz, "
                 "IsolationWidth, CollisionEnergy FROM DiaFrameMsMsWindows;"))
    {
      qDebug();
      throw ExceptionNotFound(QObject::tr("ERROR : executing %1 :\n%2")
                                .arg(query.lastQuery())
                                .arg(query.lastError().text()));
    }
  else
    {
 
      qDebug() << query.size();
      if(query.size() == 0)
        {
          qDebug();
          // FIXME: this exception is thrown, but my file is DDA.
          // That's not a nice info...
          throw ExceptionNotFound(
            QObject::tr("ERROR : DiaFrameMsMsWindows table is empty in SqlLite database"));
        }
      MsMsWindow msms_window;
      while(query.next())
        {
          qDebug();
          msms_window.WindowGroup     = query.value(0).toLongLong();
          msms_window.ScanNumBegin    = query.value(1).toLongLong();
          msms_window.ScanNumEnd      = query.value(2).toLongLong();
          msms_window.IsolationMz     = query.value(3).toDouble();
          msms_window.IsolationWidth  = query.value(4).toDouble();
          msms_window.CollisionEnergy = query.value(5).toDouble();
 
          qDebug();
          m_msMsWindowGroupList.addInGroup(msms_window);
        }
 
      if(m_msMsWindowGroupList.size() == 0)
        {
          qDebug();
          // FIXME: this exception is thrown, but my file is DDA.
          // That's not a nice info...
          throw ExceptionNotFound(
            QObject::tr("ERROR : no DiaFrameMsMsWindows in SqlLite database"));
        }
    }
 
  qDebug();
  std::size_t max_frame_id = 0;
  if(!query.exec("SELECT Frame, WindowGroup FROM DiaFrameMsMsInfo;"))
    {
      qDebug();
      throw pappso::ExceptionNotFound(
        QObject::tr("ERROR : no DiaFrameMsMsInfo in SqlLite database"));
    }
  else
    {
      qDebug() << query.size();
      if(query.size() == 0)
        {
          qDebug();
          throw pappso::ExceptionNotFound(
            QObject::tr("ERROR : no DiaFrameMsMsInfo in SqlLite database"));
        }
      while(query.next())
        {
          qDebug();
          max_frame_id = query.value(0).toLongLong();
          m_mapFrame2WindowGroupPtr[query.value(0).toLongLong()] =
            m_msMsWindowGroupList.getWindowGroupPtrByGroupId(query.value(1).toLongLong());
        }
 
      if(m_mapFrame2WindowGroupPtr.size() == 0)
        {
          qDebug();
          throw ExceptionNotFound(QObject::tr("ERROR : no DiaFrameMsMsInfo in SqlLite database"));
        }
    }
  qDebug();
 
  m_frameId2GlobalSliceIndexBegin.resize(max_frame_id + 1);
 
  m_ms2frameId2Ms1FrameId.resize(max_frame_id + 1);
  std::size_t last_ms1       = 0;
  std::size_t previous_frame = 0;
 
  std::size_t cumul_global_slice = 0;
  m_ms1frameIdList.clear();
  for(const std::pair<std::size_t, MsMsWindowGroup *> pair_frame_pgroup : m_mapFrame2WindowGroupPtr)
    {
      if(previous_frame != pair_frame_pgroup.first - 1)
        {
          last_ms1 = previous_frame + 1;
          m_ms1frameIdList.push_back(last_ms1);
        }
      previous_frame = pair_frame_pgroup.first;
      qDebug() << pair_frame_pgroup.first << " " << m_frameId2GlobalSliceIndexBegin.size();
      FrameSliceRange slice_range;
      qDebug() << pair_frame_pgroup.first;
      slice_range.frame_id                                     = pair_frame_pgroup.first;
      slice_range.begin_global_slice_index                     = cumul_global_slice;
      m_frameId2GlobalSliceIndexBegin[pair_frame_pgroup.first] = cumul_global_slice;
      m_ms2frameId2Ms1FrameId[pair_frame_pgroup.first]         = last_ms1;
 
      cumul_global_slice += pair_frame_pgroup.second->size();
 
      if(cumul_global_slice > 0)
        slice_range.end_global_slice_index = cumul_global_slice - 1;
      m_frameSliceRangeList.push_back(slice_range);
    }
 
  qDebug();
  m_totalSlicesCount = cumul_global_slice;
}
 
TimsDiaSlices::~TimsDiaSlices()
{
}
 
const std::map<std::size_t, TimsDiaSlices::MsMsWindowGroup *> &
TimsDiaSlices::getMapFrame2WindowGroupPtr() const
{
  return m_mapFrame2WindowGroupPtr;
}
const TimsDiaSlices::MsMsWindowGroupList &
TimsDiaSlices::getMsMsWindowGroupList() const
{
  return m_msMsWindowGroupList;
}
 
std::size_t
TimsDiaSlices::getGlobalSliceIndexBeginByFrameId(std::size_t frame_id) const
{
  return (m_frameId2GlobalSliceIndexBegin[frame_id]);
}
std::size_t
TimsDiaSlices::getFrameIdByGlobalSliceIndex(std::size_t global_slice_id) const
{
  return (getFrameSliceRangeByGlobalSliceIndex(global_slice_id).frame_id);
}
 
const TimsDiaSlices::FrameSliceRange &
TimsDiaSlices::getFrameSliceRangeByGlobalSliceIndex(std::size_t global_slice_index) const
{
 
  qDebug();
  auto it = std::find_if(m_frameSliceRangeList.cbegin(),
                         m_frameSliceRangeList.cend(),
                         [global_slice_index](const FrameSliceRange &x) {
                           if(global_slice_index < x.begin_global_slice_index)
                             return false;
                           if(global_slice_index > x.end_global_slice_index)
                             return false;
                           return true;
                         });
 
  if(it != m_frameSliceRangeList.cend())
    {
      return *it;
    }
  qDebug();
  throw ExceptionNotFound(
    QObject::tr("ERROR : GlobalSliceIndex %1 not found").arg(global_slice_index));
}
 
std::size_t
TimsDiaSlices::getTotalSlicesCount() const
{
  return m_totalSlicesCount;
}
 
const TimsDiaSlices::MsMsWindow &
TimsDiaSlices::getMsMsWindowByGlobalSliceIndex(std::size_t global_slice_index) const
{
  qDebug() << global_slice_index;
  const TimsDiaSlices::FrameSliceRange &range =
    getFrameSliceRangeByGlobalSliceIndex(global_slice_index);
 
  qDebug() << "range.frame_id=" << range.frame_id
           << " range.begin_global_slice_index=" << range.begin_global_slice_index;
  return m_mapFrame2WindowGroupPtr.at(range.frame_id)
    ->at(global_slice_index - range.begin_global_slice_index);
}
 
std::size_t
TimsDiaSlices::getLastMs1FrameIdByMs2FrameId(std::size_t frame_id) const
{
 
  return m_ms2frameId2Ms1FrameId.at(frame_id);
}
 
 
TimsDataFastMap &
TimsDiaSlices::getCombinedMs2ScansByGlobalSliceIndex(std::size_t global_slice_index) const
{
 
  qDebug();
  TimsDataFastMap &raw_spectrum = TimsDataFastMap::getTimsDataFastMapInstance();
  raw_spectrum.clear();
  try
    {
      FrameSliceRange range = getFrameSliceRangeByGlobalSliceIndex(global_slice_index);
 
      MsMsWindow window = m_mapFrame2WindowGroupPtr.at(range.frame_id)
                            ->at(global_slice_index - range.begin_global_slice_index);
      qDebug();
      TimsFrameCstSPtr tims_frame;
      tims_frame = mp_timsDataOrigin->getTimsFrameCstSPtrCached(range.frame_id);
      qDebug();
      /*combiner.combine(combiner_result,
        tims_frame.get()->cumulateScanToTrace(
        scan_mobility_start, scan_mobility_end));*/
      tims_frame.get()->combineScansInTofIndexIntensityMap(
        raw_spectrum, window.ScanNumBegin, window.ScanNumEnd);
      qDebug();
    }
 
  catch(PappsoException &error)
    {
      throw PappsoException(QObject::tr("ERROR in %1 (global_slice_index=%2):\n%3")
                              .arg(__FUNCTION__)
                              .arg(global_slice_index)
                              .arg(error.qwhat()));
    }
  catch(std::exception &error)
    {
      qDebug() << QString("Failure %1 ").arg(error.what());
    }
  return raw_spectrum;
  qDebug();
}
void
TimsDiaSlices::getMs2QualifiedSpectrumByGlobalSliceIndex(const MsRunIdCstSPtr &msrun_id,
                                                         QualifiedMassSpectrum &mass_spectrum,
                                                         std::size_t global_slice_index,
                                                         bool want_binary_data) const
{
 
  qDebug();
  try
    {
      FrameSliceRange range = getFrameSliceRangeByGlobalSliceIndex(global_slice_index);
 
      MsMsWindow window = m_mapFrame2WindowGroupPtr.at(range.frame_id)
                            ->at(global_slice_index - range.begin_global_slice_index);
 
      MassSpectrumId spectrum_id;
 
      spectrum_id.setSpectrumIndex(global_slice_index);
      spectrum_id.setNativeId(
        QString("global_slice_index=%1 frame=%2 begin=%3 end=%4 group=%5 slice=%6")
          .arg(global_slice_index)
          .arg(range.frame_id)
          .arg(window.ScanNumBegin)
          .arg(window.ScanNumEnd)
          .arg(window.WindowGroup)
          .arg(window.SliceIndex));
 
      spectrum_id.setMsRunId(msrun_id);
 
      mass_spectrum.setMassSpectrumId(spectrum_id);
 
      mass_spectrum.setMsLevel(2);
      // mass_spectrum.setPrecursorSpectrumIndex(spectrum_descr.ms1_index);
 
      mass_spectrum.setEmptyMassSpectrum(true);
 
      qDebug();
 
      // get precursor ion data
 
      // mass_spectrum.appendPrecursorIonData(spectrum_descr.precursor_ion_data);
 
      mass_spectrum.setPrecursorNativeId(QString("frame_id=%1 begin=%2 end=%3 group=%4 slice=%5")
                                           .arg(m_ms2frameId2Ms1FrameId.at(range.frame_id))
                                           .arg(window.ScanNumBegin)
                                           .arg(window.ScanNumEnd)
                                           .arg(window.WindowGroup)
                                           .arg(window.SliceIndex));
 
      mass_spectrum.setParameterValue(QualifiedMassSpectrumParameter::IsolationWindowTargetMz,
                                      window.IsolationMz);
      mass_spectrum.setParameterValue(QualifiedMassSpectrumParameter::IsolationWindowLowerOffset,
                                      window.IsolationWidth);
      mass_spectrum.setParameterValue(QualifiedMassSpectrumParameter::IsolationWindowUpperOffset,
                                      window.IsolationWidth);
 
      mass_spectrum.setParameterValue(QualifiedMassSpectrumParameter::CollisionEnergy,
                                      window.CollisionEnergy);
 
 
      qDebug() << window.ScanNumBegin;
      mass_spectrum.setParameterValue(QualifiedMassSpectrumParameter::TimsFrameIonMobScanIndexBegin,
                                      (qulonglong)window.ScanNumBegin);
 
 
      qDebug() << window.ScanNumEnd;
      mass_spectrum.setParameterValue(QualifiedMassSpectrumParameter::TimsFrameIonMobScanIndexEnd,
                                      (qulonglong)window.ScanNumEnd);
      qDebug();
      TimsFrameBaseCstSPtr tims_frame;
 
 
      if(want_binary_data)
        {
          qDebug() << "bindec";
          tims_frame = mp_timsDataOrigin->getTimsFrameCstSPtrCached(range.frame_id);
        }
      else
        {
          tims_frame = mp_timsDataOrigin->getTimsFrameBaseCstSPtrCached(range.frame_id);
        }
 
      mass_spectrum.setRtInSeconds(tims_frame.get()->getRtInSeconds());
 
 
      mass_spectrum.setParameterValue(
        QualifiedMassSpectrumParameter::IonMobOneOverK0Begin,
        tims_frame.get()->getOneOverK0Transformation(window.ScanNumBegin));
 
      mass_spectrum.setParameterValue(
        QualifiedMassSpectrumParameter::IonMobOneOverK0End,
        tims_frame.get()->getOneOverK0Transformation(window.ScanNumEnd));
 
      qDebug();
      /*combiner.combine(combiner_result,
        tims_frame.get()->cumulateScanToTrace(
        scan_mobility_start, scan_mobility_end));*/
      if(want_binary_data)
        {
 
          TimsDataFastMap &raw_spectrum = TimsDataFastMap::getTimsDataFastMapInstance();
          raw_spectrum.clear();
          tims_frame.get()->combineScansInTofIndexIntensityMap(
            raw_spectrum, window.ScanNumBegin, window.ScanNumEnd);
 
 
          Trace trace;
          if(m_builtinMs2Centroid)
            {
              // raw_spectrum.removeArtefactualSpike();
              raw_spectrum.builtInCentroid();
            }
 
          trace = tims_frame.get()->getTraceFromTofIndexIntensityMap(raw_spectrum);
 
          if(trace.size() > 0)
            {
              mass_spectrum.setMassSpectrumSPtr(MassSpectrum(trace).makeMassSpectrumSPtr());
              mass_spectrum.setEmptyMassSpectrum(false);
            }
          else
            {
              mass_spectrum.setMassSpectrumSPtr(nullptr);
              mass_spectrum.setEmptyMassSpectrum(true);
            }
        }
      qDebug();
    }
 
  catch(PappsoException &error)
    {
      throw PappsoException(QObject::tr("ERROR in %1 (global_slice_index=%2):\n%3")
                              .arg(__FUNCTION__)
                              .arg(global_slice_index)
                              .arg(error.qwhat()));
    }
  catch(std::exception &error)
    {
      qDebug() << QString("Failure %1 ").arg(error.what());
    }
  qDebug();
}
 
void
TimsDiaSlices::getMs1QualifiedSpectrumByGlobalSliceIndex(const MsRunIdCstSPtr &msrun_id,
                                                         QualifiedMassSpectrum &mass_spectrum,
                                                         std::size_t global_slice_index,
                                                         bool want_binary_data,
                                                         int rt_position) const
{
 
  qDebug();
  try
    {
      FrameSliceRange range = getFrameSliceRangeByGlobalSliceIndex(global_slice_index);
 
      MsMsWindow window = m_mapFrame2WindowGroupPtr.at(range.frame_id)
                            ->at(global_slice_index - range.begin_global_slice_index);
 
      std::size_t ms1_frame_id = getLastMs1FrameIdByMs2FrameId(range.frame_id);
      if(rt_position != 0)
        {
          std::size_t i = 0;
          for(i = 0; i < m_ms1frameIdList.size(); i++)
            {
              if(m_ms1frameIdList.at(i) == ms1_frame_id)
                break;
            }
          long index = (long)i + rt_position;
          if((index > -1) && (index < (long)m_ms1frameIdList.size()))
            ms1_frame_id = m_ms1frameIdList.at(index);
        }
 
      MassSpectrumId spectrum_id;
 
      spectrum_id.setSpectrumIndex(global_slice_index);
      spectrum_id.setNativeId(
        QString("global_slice_index=%1 frame=%2 begin=%3 end=%4 group=%5 slice=%6")
          .arg(global_slice_index)
          .arg(ms1_frame_id)
          .arg(window.ScanNumBegin)
          .arg(window.ScanNumEnd)
          .arg(window.WindowGroup)
          .arg(window.SliceIndex));
 
      spectrum_id.setMsRunId(msrun_id);
 
      mass_spectrum.setMassSpectrumId(spectrum_id);
 
      mass_spectrum.setMsLevel(1);
      // mass_spectrum.setPrecursorSpectrumIndex(spectrum_descr.ms1_index);
 
      mass_spectrum.setEmptyMassSpectrum(true);
 
      qDebug();
 
      // get precursor ion data
 
      // mass_spectrum.appendPrecursorIonData(spectrum_descr.precursor_ion_data);
 
 
      qDebug() << window.ScanNumBegin;
      mass_spectrum.setParameterValue(QualifiedMassSpectrumParameter::TimsFrameIonMobScanIndexBegin,
                                      (qulonglong)window.ScanNumBegin);
 
 
      qDebug() << window.ScanNumEnd;
      mass_spectrum.setParameterValue(QualifiedMassSpectrumParameter::TimsFrameIonMobScanIndexEnd,
                                      (qulonglong)window.ScanNumEnd);
      qDebug();
      TimsFrameBaseCstSPtr tims_frame;
 
 
      if(want_binary_data)
        {
          qDebug() << "bindec";
          tims_frame = mp_timsDataOrigin->getTimsFrameCstSPtrCached(ms1_frame_id);
        }
      else
        {
          tims_frame = mp_timsDataOrigin->getTimsFrameBaseCstSPtrCached(ms1_frame_id);
        }
      mass_spectrum.setRtInSeconds(tims_frame.get()->getRtInSeconds());
 
 
      mass_spectrum.setParameterValue(
        QualifiedMassSpectrumParameter::IonMobOneOverK0Begin,
        tims_frame.get()->getOneOverK0Transformation(window.ScanNumBegin));
 
      mass_spectrum.setParameterValue(
        QualifiedMassSpectrumParameter::IonMobOneOverK0End,
        tims_frame.get()->getOneOverK0Transformation(window.ScanNumEnd));
 
      qDebug();
      /*combiner.combine(combiner_result,
        tims_frame.get()->cumulateScanToTrace(
        scan_mobility_start, scan_mobility_end));*/
      if(want_binary_data)
        {
 
          double mz = (window.IsolationMz - (window.IsolationWidth / 2));
          qDebug() << mz;
 
          quint32 tof_index_begin =
            tims_frame.get()->getMzCalibrationInterfaceSPtr().get()->getTofIndexFromMz(mz);
 
          mz = (window.IsolationMz + (window.IsolationWidth / 2));
          qDebug() << mz;
          quint32 tof_index_end =
            tims_frame.get()->getMzCalibrationInterfaceSPtr().get()->getTofIndexFromMz(mz);
 
          qDebug() << "tof_index_begin=" << tof_index_begin << " tof_index_end=" << tof_index_end;
 
          TimsDataFastMap &raw_spectrum = TimsDataFastMap::getTimsDataFastMapInstance();
          raw_spectrum.clear();
          tims_frame.get()->combineScansInTofIndexIntensityMap(
            raw_spectrum, window.ScanNumBegin, window.ScanNumEnd, tof_index_begin, tof_index_end);
 
 
          Trace trace;
          if(m_builtinMs2Centroid)
            {
              // raw_spectrum.removeArtefactualSpike();
              raw_spectrum.builtInCentroid();
            }
 
          trace = tims_frame.get()->getTraceFromTofIndexIntensityMap(raw_spectrum);
 
          if(trace.size() > 0)
            {
              mass_spectrum.setMassSpectrumSPtr(MassSpectrum(trace).makeMassSpectrumSPtr());
              mass_spectrum.setEmptyMassSpectrum(false);
            }
          else
            {
              mass_spectrum.setMassSpectrumSPtr(nullptr);
              mass_spectrum.setEmptyMassSpectrum(true);
            }
        }
      qDebug();
    }
 
  catch(PappsoException &error)
    {
      throw PappsoException(QObject::tr("ERROR in %1 (global_slice_index=%2):\n%3")
                              .arg(__FUNCTION__)
                              .arg(global_slice_index)
                              .arg(error.qwhat()));
    }
  catch(std::exception &error)
    {
      qDebug() << QString("Failure %1 ").arg(error.what());
    }
  qDebug();
}
 
} // namespace pappso