mzcborprecursorfinder.cpp

Go to the documentation of this file.

/**
 * \file pappsomspp/core/processing/cbor/mzcbor/copyreplace/mzcborprecursorfinder.cpp
 * \date 2/10/2026
 * \author Olivier Langella
 * \brief read mzcbor stream or file, process it and write processed mzcbor output
 */
 
/*******************************************************************************
 * Copyright (c) 2026 Olivier Langella <Olivier.Langella@cnrs.fr>.
 *
 * This file is part of PAPPSOms-tools.
 *
 *     PAPPSOms-tools 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-tools 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-tools.  If not, see <http://www.gnu.org/licenses/>.
 *
 ******************************************************************************/
 
#include "filtercomplement.h"
#include "filterms1windowselection.h"
#include "mzcborprecursorfinder.h"
#include "pappsomspp/core/exception/exceptionnotfound.h"
#include "pappsomspp/core/mzrange.h"
#include "pappsomspp/core/precision.h"
#include "pappsomspp/core/processing/cbor/mzcbor/copyreplace/mzcborcopyreplacebase.h"
#include "pappsomspp/core/pappsoexception.h"
#include "pappsomspp/core/processing/cbor/mzcbor/copyreplace/qcborspectrum.h"
#include "pappsomspp/core/processing/cbor/mzcbor/copyreplace/qcborspectrumprecursor.h"
#include "pappsomspp/core/processing/cbor/mzcbor/precursor.h"
#include "pappsomspp/core/processing/filters/filterpass.h"
#include "pappsomspp/core/trace/trace.h"
#include <memory>
#include <pappsomspp/core/processing/filters/filterchargedeconvolution.h>
#include <pappsomspp/core/processing/filters/filterresample.h>
#include <algorithm>
#include <qjsonarray.h>
#include <qlogging.h>
#include <qobject.h>
namespace pappso
{
namespace cbor
{
namespace mzcbor
{
MzcborPrecursorFinder::MzcborPrecursorFinder(CborStreamWriter *cbor_output_p,
                                             const QJsonObject &json_params)
  : MzcborCopyReplaceBase(cbor_output_p, "precursor_finder")
{
  m_parameters = json_params;
 
 
  if(getParameterValue("ms1_tolerance_unit").toString() == "dalton")
    {
      m_ms1Tolerance =
        pappso::PrecisionFactory::getDaltonInstance(getParameterValue("ms1_tolerance").toDouble());
    }
  else if(getParameterValue("ms1_tolerance_unit").toString() == "ppm")
    {
      m_ms1Tolerance =
        pappso::PrecisionFactory::getPpmInstance(getParameterValue("ms1_tolerance").toDouble());
    }
  else
    {
      throw pappso::ExceptionNotFound(QObject::tr("ms1_tolerance_unit %1 not known")
                                        .arg(getParameterValue("ms1_tolerance_unit").toString()));
    }
 
  if(getParameterValue("ms2_tolerance_unit").toString() == "dalton")
    {
      m_ms2Tolerance =
        pappso::PrecisionFactory::getDaltonInstance(getParameterValue("ms2_tolerance").toDouble());
    }
  else if(getParameterValue("ms2_tolerance_unit").toString() == "ppm")
    {
      m_ms2Tolerance =
        pappso::PrecisionFactory::getPpmInstance(getParameterValue("ms2_tolerance").toDouble());
    }
  else
    {
      throw pappso::ExceptionNotFound(QObject::tr("ms2_tolerance_unit %1 not known")
                                        .arg(getParameterValue("ms2_tolerance_unit").toString()));
    }
 
  m_minimumComplementIon    = getParameterValue("minimum_complement_ion").toInt();
  m_intensityRatioThreshold = getParameterValue("intensity_ratio_threshold").toDouble();
 
  m_ms2nMostIntense = getParameterValue("ms2_max_peaks").toInt();
 
  m_potentialPrecursorChargeStateList.clear();
  QJsonArray charge_arr = getParameterValue("potential_charge_states").toArray();
  for(auto charge : charge_arr)
    {
      m_potentialPrecursorChargeStateList.push_back(charge.toInt());
    }
}
MzcborPrecursorFinder::~MzcborPrecursorFinder()
{
}
 
 
QJsonValue
MzcborPrecursorFinder::getParameterValue(const QString &name) const
{
  QJsonValue value;
  if(m_parameters.contains(name))
    {
      return m_parameters.value(name);
    }
  else
    {
      throw pappso::ExceptionNotFound(QObject::tr("parameter %1 not defined precursor finder JSON "
                                                  "parameters")
                                        .arg(name));
    }
  return value;
}
 
void
MzcborPrecursorFinder::setRetentionTimeRangeInSeconds(double rt_start, double rt_stop)
{
  m_rtStart = rt_start;
  m_rtStop  = rt_stop;
}
void
MzcborPrecursorFinder::processQcborSpectrum(SpectrumTube &spectrum_tube,
                                            QcborSpectrum &cbor_spectrum)
{
  if(cbor_spectrum.getMsLevel() == 1)
    {
      QString old_ms1_id = cbor_spectrum.getId();
      m_spectrumMs1List.push_back({old_ms1_id, cbor_spectrum});
      if(m_spectrumMs1List.size() > m_maxSpectrumMs1ListSize)
        {
          m_spectrumMs1List.pop_front();
        }
      QString new_ms1_id = spectrum_tube.push_back(cbor_spectrum);
 
      m_ms1idOldToNew.insert({old_ms1_id, new_ms1_id});
    }
  else
    {
      double rt = cbor_spectrum.getRetentionTimeInSeconds();
      if((rt > m_rtStart) && (rt < m_rtStop))
        {
 
          std::vector<QcborSpectrumPrecursor> precursor_list = cbor_spectrum.getPrecursorList();
          if(precursor_list.size() == 1)
            {
              //"spectrumRef":"controllerType=0 controllerNumber=1 scan=2711",
              QString spectrum_ref = precursor_list.at(0).getSpectrumRef();
 
              auto it_ms1 =
                std::find_if(m_spectrumMs1List.begin(),
                             m_spectrumMs1List.end(),
                             [spectrum_ref](const std::pair<QString, QcborSpectrum> &pair) {
                               return pair.first == spectrum_ref;
                             });
              if(it_ms1 != m_spectrumMs1List.end())
                {
 
                  // replace with new spectrum_ref
                  QString new_spectrum_ref = m_ms1idOldToNew.at(spectrum_ref);
 
 
                  QcborSpectrumPrecursor precursor = precursor_list.at(0);
 
                  pappso::MzRange original_ms1_peak_range(
                    precursor.getSelectedIonList().at(0).getMz(), m_ms1Tolerance);
                  // qWarning() << "original precursor charge="
                  //            << precursor.getSelectedIonList().at(0).getChargeState();
                  precursor.setSpectrumRef(new_spectrum_ref);
                  std::vector<QcborSpectrumPrecursor> new_precursor_list;
                  new_precursor_list.push_back(precursor);
                  cbor_spectrum.setPrecursorList(new_precursor_list);
 
 
                  std::vector<QcborSpectrum> spectrum_list =
                    findPotentialPrecursorIonData(it_ms1->second, cbor_spectrum, precursor);
 
                  int i                         = 0;
                  bool report_original_spectrum = true;
                  for(QcborSpectrum &spectrum : spectrum_list)
                    {
                      QcborSpectrumPrecursor candidate_precursor =
                        spectrum.getPrecursorList().at(0);
 
                      // qWarning() << "candidate_precursor charge="
                      //            <<
                      //            candidate_precursor.getSelectedIonList().at(0).getChargeState();
                      if(candidate_precursor.getSelectedIonList().at(0).getChargeState() ==
                         precursor.getSelectedIonList().at(0).getChargeState())
                        {
                          if(original_ms1_peak_range.contains(
                               candidate_precursor.getSelectedIonList().at(0).getMz()))
                            {
                              report_original_spectrum = false;
                            }
                        }
 
                      QString new_id = spectrum_tube.push_back(spectrum);
                      // qWarning() << " id=" << new_id << " candidate_precursor charge="
                      //            << spectrum.getPrecursorList()
                      //                 .at(0)
                      //                 .getSelectedIonList()
                      //                 .at(0)
                      //                 .getChargeState();
                      i++;
                    }
                  if(report_original_spectrum)
                    {
                      QString new_id = spectrum_tube.push_back(cbor_spectrum);
                      // qWarning() << "orig id=" << new_id
                      //            << " spectrum_list.size()=" << spectrum_list.size()
                      //            << " candidate_precursor charge="
                      //            << cbor_spectrum.getPrecursorList()
                      //                 .at(0)
                      //                 .getSelectedIonList()
                      //                 .at(0)
                      //                 .getChargeState();
                    }
                }
              else
                {
                  throw pappso::ExceptionNotFound(
                    QObject::tr("ERROR: precursor ref %1 not found for current spectrum %2")
                      .arg(spectrum_ref)
                      .arg(cbor_spectrum.getId()));
                }
            }
          else
            {
              qDebug() << cbor_spectrum.getId();
              throw pappso::PappsoException(
                QObject::tr("ERROR: precursor_list.size() != 1 (%1)").arg(precursor_list.size()));
            }
        }
    }
}
 
 
std::vector<QcborSpectrum>
MzcborPrecursorFinder::findPotentialPrecursorIonData(const QcborSpectrum &ms1_spectrum,
                                                     const QcborSpectrum &ms2_spectrum,
                                                     const QcborSpectrumPrecursor &precursor)
{
 
  pappso::TraceSPtr ms1_trace_sp;
 
  auto it_cache_ms1_trace =
    std::find_if(m_cacheMs1Trace.cbegin(),
                 m_cacheMs1Trace.cend(),
                 [ms1_spectrum](const auto &item) { return item.first == &ms1_spectrum; });
 
  if(it_cache_ms1_trace == m_cacheMs1Trace.end())
    {
      ms1_trace_sp = std::make_shared<Trace>();
      ms1_spectrum.decodeTrace(*ms1_trace_sp.get());
      m_cacheMs1Trace.push_back({&ms1_spectrum, ms1_trace_sp});
      if(m_cacheMs1Trace.size() > m_maxSpectrumMs1ListSize)
        {
          m_cacheMs1Trace.pop_front();
        }
    }
  else
    {
      ms1_trace_sp = it_cache_ms1_trace->second;
    }
  pappso::Trace ms1_trace(*ms1_trace_sp.get());
  pappso::Trace ms2_trace;
  ms2_spectrum.decodeTrace(ms2_trace);
 
  // MS2 filters
  pappso::FilterChargeDeconvolution(m_ms2Tolerance).filter(ms2_trace);
  pappso::FilterResampleKeepGreater(m_ms2minimumMz).filter(ms2_trace);
  pappso::FilterGreatestY(m_ms2nMostIntense).filter(ms2_trace);
  // END MS2 filters
 
 
  std::vector<QcborSpectrum> potential_precursors;
  double precursor_intensity                 = 0;
  std::vector<SelectedIon> selected_ion_list = precursor.getSelectedIonList();
  if(selected_ion_list.size() != 1)
    {
      throw pappso::PappsoException("ERROR: selected_ion_list.size() != 1");
    }
  bool is_ok          = false;
  precursor_intensity = selected_ion_list.at(0).getIntensity(&is_ok);
  if(!is_ok)
    {
      throw pappso::PappsoException(
        QObject::tr("ERROR: precursor intensity not found in spectrum %1")
          .arg(ms2_spectrum.getId()));
      precursor_intensity = 0;
    }
 
  qDebug();
  FilterMs1WindowSelection window(precursor);
  qDebug();
  window.filter(ms1_trace);
  pappso::FilterHighPass(precursor_intensity * m_intensityRatioThreshold).filter(ms1_trace);
 
  for(const pappso::DataPoint &data_point : ms1_trace)
    {
      for(int charge : m_potentialPrecursorChargeStateList)
        {
          double intensity = data_point.y;
          double mz        = data_point.x;
 
          FilterComplement complement(mz * (double)charge, m_ms2Tolerance);
          pappso::Trace ms2(ms2_trace);
          complement.filter(ms2);
          if(ms2.size() > (std::size_t)m_minimumComplementIon)
            {
              QcborSpectrum new_ms2(ms2_spectrum);
 
              QcborSpectrumPrecursor copy_precursor(precursor);
              copy_precursor.setSelectedIon(mz, intensity, charge);
 
              std::vector<QcborSpectrumPrecursor> precursor_list;
              precursor_list.push_back(copy_precursor);
 
              qDebug() << copy_precursor;
              new_ms2.setPrecursorList(precursor_list);
 
              qDebug() << new_ms2;
              potential_precursors.push_back(new_ms2);
            }
        }
    }
 
 
  return potential_precursors;
}
 
} // namespace mzcbor
} // namespace cbor
} // namespace pappso