msrunretentiontime.cpp

Go to the documentation of this file.

 
/*******************************************************************************
 * 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/>.
 *
 * Contributors:
 *     Olivier Langella <Olivier.Langella@u-psud.fr> - initial API and
 *implementation
 ******************************************************************************/
 
#include "msrunretentiontime.h"
#include "pappsomspp/core/exception/exceptionnotpossible.h"
#include <QDebug>
#include <QObject>
 
// #include <odsstream/odsexception.h>
// #include <odsstream/odsdocwriter.h>
// #include <QFileInfo>
 
using namespace pappso;
 
template <class T>
MsRunRetentionTime<T>::MsRunRetentionTime(const std::vector<double> &msrun_retention_time_line)
  : m_ms2MedianFilter(10), m_ms2MeanFilter(15), m_ms1MeanFilter(1)
{
  m_ms1RetentionTimeVector = msrun_retention_time_line;
 
 
  std::sort(m_ms1RetentionTimeVector.begin(),
            m_ms1RetentionTimeVector.end(),
            [](const double &a, const double &b) { return (a < b); });
}
 
template <class T>
MsRunRetentionTime<T>::MsRunRetentionTime(const MsRunRetentionTime<T> &other)
  : m_ms2MedianFilter(other.m_ms2MedianFilter),
    m_ms2MeanFilter(other.m_ms2MeanFilter),
    m_ms1MeanFilter(other.m_ms1MeanFilter)
{
  m_ms1RetentionTimeVector     = other.m_ms1RetentionTimeVector;
  m_alignedRetentionTimeVector = other.m_alignedRetentionTimeVector;
 
  m_seamarks        = other.m_seamarks;
  m_valuesCorrected = other.m_valuesCorrected;
 
  m_allMs2Points = other.m_allMs2Points;
 
  m_retentionTimeReferenceMethod = other.m_retentionTimeReferenceMethod;
}
 
template <class T>
MsRunRetentionTime<T>::~MsRunRetentionTime()
{
}
 
template <typename T>
const pappso::FilterMorphoMedian &
pappso::MsRunRetentionTime<T>::getMs2MedianFilter() const
{
  return m_ms2MedianFilter;
}
 
template <class T>
void
MsRunRetentionTime<T>::setMs2MedianFilter(const FilterMorphoMedian &ms2MedianFilter)
{
  m_ms2MedianFilter = ms2MedianFilter;
}
 
template <typename T>
const pappso::FilterMorphoMean &
pappso::MsRunRetentionTime<T>::getMs2MeanFilter() const
{
  return m_ms2MeanFilter;
}
 
 
template <class T>
void
MsRunRetentionTime<T>::setMs2MeanFilter(const FilterMorphoMean &ms2MeanFilter)
{
  m_ms2MeanFilter = ms2MeanFilter;
}
 
template <typename T>
const pappso::FilterMorphoMean &
pappso::MsRunRetentionTime<T>::getMs1MeanFilter() const
{
  return m_ms1MeanFilter;
}
 
template <class T>
void
MsRunRetentionTime<T>::setMs1MeanFilter(const FilterMorphoMean &ms1MeanFilter)
{
  m_ms1MeanFilter = ms1MeanFilter;
}
 
template <class T>
const std::vector<MsRunRetentionTimeSeamarkPoint<T>> &
MsRunRetentionTime<T>::getSeamarks() const
{
  qDebug();
  return m_seamarks;
}
 
template <class T>
const std::vector<double> &
MsRunRetentionTime<T>::getAlignedRetentionTimeVector() const
{
  return m_alignedRetentionTimeVector;
}
 
template <class T>
std::size_t
MsRunRetentionTime<T>::getNumberOfCorrectedValues() const
{
  return m_valuesCorrected;
}
template <class T>
const std::vector<double> &
MsRunRetentionTime<T>::getMs1RetentionTimeVector() const
{
  return m_ms1RetentionTimeVector;
}
 
template <class T>
Trace
MsRunRetentionTime<T>::getCommonDeltaRt(
  const std::vector<MsRunRetentionTimeSeamarkPoint<T>> &other_seamarks) const
{
  Trace common_points;
  getCommonDeltaRt(common_points, other_seamarks);
  return common_points;
}
 
template <class T>
void
MsRunRetentionTime<T>::addPeptideAsSeamark(const T &peptide_id,
                                           double retentionTime,
                                           double precursorIntensity)
{
  PeptideMs2Point ms2point;
  ms2point.entityHash         = peptide_id;
  ms2point.precursorIntensity = precursorIntensity;
  ms2point.retentionTime      = retentionTime;
  m_allMs2Points.push_back(ms2point);
}
 
 
template <class T>
void
MsRunRetentionTime<T>::computeSeamarks()
{
 
  qDebug();
  if(m_allMs2Points.size() == 0)
    {
      // already computed
      return;
    }
  m_seamarks.clear();
  if(m_retentionTimeReferenceMethod == ComputeRetentionTimeReference::maximum_intensity)
    {
 
 
      std::sort(m_allMs2Points.begin(),
                m_allMs2Points.end(),
                [](const PeptideMs2Point &a, const PeptideMs2Point &b) {
                  if(a.entityHash == b.entityHash)
                    {
                      return (a.precursorIntensity > b.precursorIntensity);
                    }
                  return (a.entityHash < b.entityHash);
                });
 
      auto itend = std::unique(m_allMs2Points.begin(),
                               m_allMs2Points.end(),
                               [](const PeptideMs2Point &a, const PeptideMs2Point &b) {
                                 return (a.entityHash == b.entityHash);
                               });
 
      auto it = m_allMs2Points.begin();
      while(it != itend)
        {
          m_seamarks.push_back({it->entityHash, it->retentionTime, it->precursorIntensity});
          it++;
        }
    }
  m_allMs2Points.clear();
 
  std::sort(
    m_seamarks.begin(),
    m_seamarks.end(),
    [](const MsRunRetentionTimeSeamarkPoint<T> &a, const MsRunRetentionTimeSeamarkPoint<T> &b) {
      return (a.entityHash < b.entityHash);
    });
  qDebug();
}
 
template <class T>
void
MsRunRetentionTime<T>::getCommonDeltaRt(
  Trace &delta_rt, const std::vector<MsRunRetentionTimeSeamarkPoint<T>> &other_seamarks) const
{
 
  qDebug();
  auto it = other_seamarks.begin();
 
  for(const MsRunRetentionTimeSeamarkPoint<T> &seamark : m_seamarks)
    {
      while((it != other_seamarks.end()) && (it->entityHash < seamark.entityHash))
        {
          it++;
        }
      if(it == other_seamarks.end())
        break;
      if(it->entityHash == seamark.entityHash)
        {
          delta_rt.push_back(
            DataPoint(seamark.retentionTime, seamark.retentionTime - it->retentionTime));
        }
    }
 
  qDebug();
  if((m_ms2MedianFilter.getHalfWindowSize() * 2 + 1) >= delta_rt.size())
    {
      throw ExceptionNotPossible(QObject::tr("ERROR : MS2 alignment not possible : "
                                             "\ntoo few MS2 points (%1) in common")
                                   .arg(delta_rt.size()));
    }
 
  qDebug();
  if((m_ms2MeanFilter.getHalfWindowSize() * 2 + 1) >= delta_rt.size())
    {
      throw ExceptionNotPossible(QObject::tr("ERROR : MS2 alignment not possible : "
                                             "\ntoo few MS2 points (%1) in common")
                                   .arg(delta_rt.size()));
    }
  delta_rt.sortX();
 
  // there can be multiple entities (peptides) at one retention time
  // in this case, avoid retention time redundancy by applying unique on trace :
  delta_rt.unique();
 
  qDebug();
}
 
 
template <class T>
double
MsRunRetentionTime<T>::getFrontRetentionTimeReference() const
{
  if(isAligned())
    {
      return m_alignedRetentionTimeVector.front();
    }
  return m_ms1RetentionTimeVector.front();
}
template <class T>
double
MsRunRetentionTime<T>::getBackRetentionTimeReference() const
{
  if(isAligned())
    {
      return m_alignedRetentionTimeVector.back();
    }
  return m_ms1RetentionTimeVector.back();
}
 
 
template <class T>
double
MsRunRetentionTime<T>::translateOriginal2AlignedRetentionTime(double original_retention_time) const
{
  if(m_alignedRetentionTimeVector.size() < 3)
    {
      throw ExceptionNotPossible(QObject::tr("ERROR : too few aligned points to compute aligned "
                                             "retention time (%1)")
                                   .arg(m_ms1RetentionTimeVector.size()));
    }
  if(m_alignedRetentionTimeVector.size() != m_ms1RetentionTimeVector.size())
    {
      throw ExceptionNotPossible(QObject::tr("ERROR : m_alignedRetentionTimeVector.size() %1 != %2 "
                                             "m_ms1RetentionTimeVector.size()")
                                   .arg(m_alignedRetentionTimeVector.size())
                                   .arg(m_ms1RetentionTimeVector.size()));
    }
  auto it_plus = std::find_if(m_ms1RetentionTimeVector.begin(),
                              m_ms1RetentionTimeVector.end(),
                              [original_retention_time](const double &rt_point) {
                                return original_retention_time < rt_point;
                              });
  double rt1_a, rt2_a, rt1_b, rt2_b;
  if(it_plus == m_ms1RetentionTimeVector.end())
    {
      it_plus--;
    }
  if(it_plus == m_ms1RetentionTimeVector.begin())
    {
      it_plus++;
    }
  auto it_minus = it_plus - 1;
 
  rt1_a = *it_minus;
  rt2_a = *it_plus;
 
  double ratio = (original_retention_time - rt1_a) / (rt2_a - rt1_a);
 
  auto itref = m_alignedRetentionTimeVector.begin() +
               std::distance(m_ms1RetentionTimeVector.begin(), it_minus);
 
  rt1_b = *itref;
  itref++;
  rt2_b = *itref;
 
  return (((rt2_b - rt1_b) * ratio) + rt1_b);
}
 
template <class T>
double
MsRunRetentionTime<T>::translateAligned2OriginalRetentionTime(double aligned_retention_time) const
{
  if(m_alignedRetentionTimeVector.size() < 3)
    {
      throw ExceptionNotPossible(QObject::tr("ERROR : too few aligned points to compute aligned "
                                             "retention time (%1)")
                                   .arg(m_ms1RetentionTimeVector.size()));
    }
  if(m_alignedRetentionTimeVector.size() != m_ms1RetentionTimeVector.size())
    {
      throw ExceptionNotPossible(QObject::tr("ERROR : m_alignedRetentionTimeVector.size() %1 != %2 "
                                             "m_ms1RetentionTimeVector.size()")
                                   .arg(m_alignedRetentionTimeVector.size())
                                   .arg(m_ms1RetentionTimeVector.size()));
    }
  auto it_plus = std::find_if(
    m_alignedRetentionTimeVector.begin(),
    m_alignedRetentionTimeVector.end(),
    [aligned_retention_time](const double &rt_point) { return aligned_retention_time < rt_point; });
  double rt1_a, rt2_a, rt1_b, rt2_b;
  if(it_plus == m_alignedRetentionTimeVector.end())
    {
      it_plus--;
    }
  if(it_plus == m_alignedRetentionTimeVector.begin())
    {
      it_plus++;
    }
  auto it_minus = it_plus - 1;
 
  rt1_a = *it_minus;
  rt2_a = *it_plus;
 
  double ratio = (aligned_retention_time - rt1_a) / (rt2_a - rt1_a);
 
  auto itref = m_ms1RetentionTimeVector.begin() +
               std::distance(m_alignedRetentionTimeVector.begin(), it_minus);
 
  rt1_b = *itref;
  itref++;
  rt2_b = *itref;
 
  return (((rt2_b - rt1_b) * ratio) + rt1_b);
}
 
template <class T>
const std::vector<MsRunRetentionTimeSeamarkPoint<T>>
MsRunRetentionTime<T>::getSeamarksReferences() const
{
  std::vector<MsRunRetentionTimeSeamarkPoint<T>> other_seamarks = m_seamarks;
  for(auto &seamark : other_seamarks)
    {
      seamark.retentionTime = translateOriginal2AlignedRetentionTime(seamark.retentionTime);
    }
  return other_seamarks;
}
 
template <class T>
bool
MsRunRetentionTime<T>::isAligned() const
{
  return (m_alignedRetentionTimeVector.size() > 0);
}
 
template <class T>
Trace
MsRunRetentionTime<T>::align(const MsRunRetentionTime<T> &msrun_retention_time_reference)
{
  computeSeamarks();
  std::vector<MsRunRetentionTimeSeamarkPoint<T>> other_seamarks;
  if(msrun_retention_time_reference.isAligned())
    {
      other_seamarks = msrun_retention_time_reference.getSeamarksReferences();
    }
  else
    {
      other_seamarks = msrun_retention_time_reference.getSeamarks();
    }
  qDebug();
  if((m_ms1MeanFilter.getHalfWindowSize() * 2 + 1) >= m_ms1RetentionTimeVector.size())
    {
      throw ExceptionNotPossible(QObject::tr("ERROR : MS1 alignment not possible : "
                                             "\ntoo few MS1 points (%1)")
                                   .arg(m_ms1RetentionTimeVector.size()));
    }
 
  qDebug() << m_seamarks[0].entityHash << " " << m_seamarks[0].retentionTime << " "
           << other_seamarks[0].entityHash << other_seamarks[0].retentionTime << " ";
  // both seamarks has to be ordered
  Trace common_points;
  getCommonDeltaRt(common_points, other_seamarks);
 
  // writeTrace("lib_ms2_delta_rt.ods", common_points);
 
  qDebug() << common_points.front().x << " " << common_points.front().y;
  m_ms2MedianFilter.filter(common_points);
  // writeTrace("lib_ms2_delta_rt_median.ods", common_points);
  m_ms2MeanFilter.filter(common_points);
  // writeTrace("lib_ms2_delta_rt_mean.ods", common_points);
  // convert common delta rt to real retention times (for convenience)
  qDebug() << common_points.front().x << " " << common_points.front().y;
 
 
  // add a first point to ensure coherence:
  DataPoint first_point;
  first_point.x = m_ms1RetentionTimeVector.front() - (double)1;
  if(first_point.x < 0)
    {
      first_point.x = 0;
    }
  first_point.y = m_ms1RetentionTimeVector.front() -
                  msrun_retention_time_reference.getFrontRetentionTimeReference();
 
  common_points.push_back(first_point);
  // add a last point to ensure coherence:
  DataPoint last_point;
  last_point.x = m_ms1RetentionTimeVector.back() + 1;
  last_point.y = m_ms1RetentionTimeVector.back() -
                 msrun_retention_time_reference.getBackRetentionTimeReference();
  common_points.push_back(last_point);
  common_points.sortX();
 
  // now, it is possible for each time range to give a new MS1 time using a
  // linear regression on MS2 corrected times
  m_alignedRetentionTimeVector.clear();
 
  qDebug() << common_points.front().x << " " << common_points.front().y;
 
  Trace ms1_aligned_points;
 
  linearRegressionMs2toMs1(ms1_aligned_points, common_points);
 
  // writeTrace("lib_ms1_map_rt.ods", ms1_aligned_points);
  qDebug();
  // smoothing on MS1 points
  m_ms1MeanFilter.filter(ms1_aligned_points);
 
  // writeTrace("lib_ms1_map_rt_mean.ods", ms1_aligned_points);
  // final aligned retentionTime vector
 
  for(DataPoint &data_point : ms1_aligned_points)
    {
      data_point.y = (data_point.x - data_point.y);
    }
 
  qDebug();
  // Here, the correction parameter is the slope of old rt points curve
  // (divided by 4 to get a finer correction).
  double correction_parameter =
    (m_ms1RetentionTimeVector.back() - m_ms1RetentionTimeVector.front()) /
    (ms1_aligned_points.size());
  // set_correction_parameter(correction_parameter / 4);
  correction_parameter = correction_parameter / (double)4;
  correctNewTimeValues(ms1_aligned_points, correction_parameter);
 
  m_alignedRetentionTimeVector = ms1_aligned_points.yValues();
 
  qDebug();
  return ms1_aligned_points;
}
 
 
template <class T>
void
MsRunRetentionTime<T>::linearRegressionMs2toMs1(Trace &ms1_aligned_points,
                                                const Trace &common_points)
{
 
  // first slope :
  std::vector<DataPoint>::const_iterator itms2     = common_points.begin();
  std::vector<DataPoint>::const_iterator itms2next = itms2 + 1;
  if(itms2next == common_points.end())
    {
      // error
      throw ExceptionNotPossible(QObject::tr("ERROR : MS1 alignment not possible : "
                                             "\ntoo few common points (%1)")
                                   .arg(common_points.size()));
    }
  qDebug() << "() itms2->x=" << itms2->x << " itms2->y=" << itms2->y;
 
  for(double &original_rt_point : m_ms1RetentionTimeVector)
    {
      DataPoint ms1_point;
      ms1_point.x = original_rt_point;
 
      while(ms1_point.x > itms2next->x)
        {
          itms2++;
          itms2next++;
        }
 
      double ratio = (itms2next->x - itms2->x);
      if(ratio != 0)
        {
          ratio = (ms1_point.x - itms2->x) / ratio;
        }
      else
        {
          // avoid division by zero
          ratio = 1;
        }
      // qDebug() << __FILE__ << "@" << __LINE__ << __FUNCTION__ << "() " <<
      // ratio;
 
      ms1_point.y = itms2->y + ((itms2next->y - itms2->y) * ratio);
 
      // qDebug() << __FILE__ << "@" << __LINE__ << __FUNCTION__ << "() "
      //         << ms1_point.y;
      ms1_aligned_points.push_back(ms1_point);
    }
}
 
template <class T>
void
MsRunRetentionTime<T>::correctNewTimeValues(Trace &ms1_aligned_points, double correction_parameter)
{
 
  m_valuesCorrected = 0;
  auto new_it(ms1_aligned_points.begin());
  auto new_nextit(ms1_aligned_points.begin());
  new_nextit++;
  for(; new_nextit != ms1_aligned_points.end(); ++new_nextit, ++new_it)
    {
      if(new_nextit->y < new_it->y)
        {
          ++m_valuesCorrected;
          new_nextit->y = new_it->y + correction_parameter;
        }
    }
}
 
template <class T>
void
pappso::MsRunRetentionTime<T>::setAlignedRetentionTimeVector(
  const std::vector<double> &aligned_times)
{
 
  if(aligned_times.size() == m_ms1RetentionTimeVector.size())
    {
      m_alignedRetentionTimeVector = aligned_times;
    }
  else
    {
      if(aligned_times.size() == m_ms1RetentionTimeVector.size() * 2)
        {
          m_alignedRetentionTimeVector = m_ms1RetentionTimeVector;
          for(std::size_t i = 0; i < m_ms1RetentionTimeVector.size(); i++)
            {
              if(aligned_times[2 * i] != m_ms1RetentionTimeVector[i])
                {
                  throw pappso::PappsoException(
                    QObject::tr("ERROR : aligned_times (size=%1) vector does not have "
                                "required size (size=%2)")
                      .arg(aligned_times.size())
                      .arg(m_ms1RetentionTimeVector.size()));
                }
              m_alignedRetentionTimeVector[i] = aligned_times[(2 * i) + 1];
            }
        }
      else
        {
          throw ExceptionNotPossible(
            QObject::tr("ERROR : aligned_times (size=%1) vector does not have "
                        "required size (size=%2)")
              .arg(aligned_times.size())
              .arg(m_ms1RetentionTimeVector.size()));
        }
    }
}
 
 
template <class T>
Trace
MsRunRetentionTime<T>::getCommonSeamarksDeltaRt(
  const MsRunRetentionTime<T> &msrun_retention_time_reference) const
{
  // computeSeamarks();
  std::vector<MsRunRetentionTimeSeamarkPoint<T>> other_seamarks;
  if(msrun_retention_time_reference.isAligned())
    {
      other_seamarks = msrun_retention_time_reference.getSeamarksReferences();
    }
  else
    {
      other_seamarks = msrun_retention_time_reference.getSeamarks();
    }
  qDebug();
 
  qDebug() << m_seamarks[0].entityHash << " " << m_seamarks[0].retentionTime << " "
           << other_seamarks[0].entityHash << other_seamarks[0].retentionTime << " ";
  // both seamarks has to be ordered
  Trace common_points;
  getCommonDeltaRt(common_points, other_seamarks);
  return common_points;
}