timsframe.cpp

Go to the documentation of this file.

/**
 * \file pappsomspp/vendors/tims/timsframe.cpp
 * \date 23/08/2019
 * \author Olivier Langella
 * \brief handle a single Bruker's TimsTof frame
 */
 
/*******************************************************************************
 * 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 "timsframe.h"
#include "../../../pappsomspp/pappsoexception.h"
#include "../../../pappsomspp/exception/exceptionoutofrange.h"
#include <QDebug>
#include <QtEndian>
#include <memory>
#include <solvers.h>
#include "xicextractor/timsdirectxicextractor.h"
 
 
namespace pappso
{
 
TimsFrame::TimsFrame(std::size_t timsId,
                     quint32 scanNum,
                     char *p_bytes,
                     std::size_t len)
  : TimsFrameBase(timsId, scanNum)
{
  // langella@themis:~/developpement/git/bruker/cbuild$
  // ./src/sample/timsdataSamplePappso
  // /gorgone/pappso/fichiers_fabricants/Bruker/Demo_TimsTOF_juin2019/Samples/1922001/1922001-1_S-415_Pep_Pur-1ul_Slot1-10_1_2088.d/
  qDebug() << timsId;
 
  m_timsDataFrame.resize(len);
 
  if(p_bytes != nullptr)
    {
      unshufflePacket(p_bytes);
    }
  else
    {
      if(m_scanNumber == 0)
        {
 
          throw pappso::PappsoException(
            QObject::tr("TimsFrame::TimsFrame(%1,%2,nullptr,%3) FAILED")
              .arg(m_timsId)
              .arg(m_scanNumber)
              .arg(len));
        }
    }
}
 
TimsFrame::TimsFrame(const TimsFrame &other) : TimsFrameBase(other)
{
}
 
TimsFrame::~TimsFrame()
{
}
 
 
void
TimsFrame::unshufflePacket(const char *src)
{
  qDebug();
  quint64 len = m_timsDataFrame.size();
  if(len % 4 != 0)
    {
      throw pappso::PappsoException(
        QObject::tr("TimsFrame::unshufflePacket error: len%4 != 0"));
    }
 
  quint64 nb_uint4 = len / 4;
 
  char *dest         = m_timsDataFrame.data();
  quint64 src_offset = 0;
 
  for(quint64 j = 0; j < 4; j++)
    {
      for(quint64 i = 0; i < nb_uint4; i++)
        {
          dest[(i * 4) + j] = src[src_offset];
          src_offset++;
        }
    }
  qDebug();
}
 
std::size_t
TimsFrame::getNbrPeaks(std::size_t scanNum) const
{
  if(m_timsDataFrame.size() == 0)
    return 0;
  /*
    if(scanNum == 0)
      {
        quint32 res = (*(quint32 *)(m_timsDataFrame.constData() + 4)) -
                      (*(quint32 *)(m_timsDataFrame.constData()-4));
        return res / 2;
      }*/
  if(scanNum == (m_scanNumber - 1))
    {
      auto nb_uint4 = m_timsDataFrame.size() / 4;
 
      std::size_t cumul = 0;
      for(quint32 i = 0; i < m_scanNumber; i++)
        {
          cumul += (*(quint32 *)(m_timsDataFrame.constData() + (i * 4)));
        }
      return (nb_uint4 - cumul) / 2;
    }
  checkScanNum(scanNum);
 
  // quint32 *res = (quint32 *)(m_timsDataFrame.constData() + (scanNum * 4));
  // qDebug() << " res=" << *res;
  return (*(quint32 *)(m_timsDataFrame.constData() + ((scanNum + 1) * 4))) / 2;
}
 
std::size_t
TimsFrame::getScanOffset(std::size_t scanNum) const
{
  std::size_t offset = 0;
  for(std::size_t i = 0; i < (scanNum + 1); i++)
    {
      offset += (*(quint32 *)(m_timsDataFrame.constData() + (i * 4)));
    }
  return offset;
}
 
 
std::vector<quint32>
TimsFrame::getScanIndexList(std::size_t scanNum) const
{
  qDebug();
  checkScanNum(scanNum);
  std::vector<quint32> scan_tof;
 
  if(m_timsDataFrame.size() == 0)
    return scan_tof;
  scan_tof.resize(getNbrPeaks(scanNum));
 
  std::size_t offset = getScanOffset(scanNum);
 
  qint32 previous = -1;
  for(std::size_t i = 0; i < scan_tof.size(); i++)
    {
      scan_tof[i] =
        (*(quint32 *)(m_timsDataFrame.constData() + (offset * 4) + (i * 8))) +
        previous;
      previous = scan_tof[i];
    }
  qDebug();
  return scan_tof;
}
 
std::vector<quint32>
TimsFrame::getScanIntensities(std::size_t scanNum) const
{
  qDebug();
  checkScanNum(scanNum);
  std::vector<quint32> scan_intensities;
 
  if(m_timsDataFrame.size() == 0)
    return scan_intensities;
 
  scan_intensities.resize(getNbrPeaks(scanNum));
 
  std::size_t offset = getScanOffset(scanNum);
 
  for(std::size_t i = 0; i < scan_intensities.size(); i++)
    {
      scan_intensities[i] = (*(quint32 *)(m_timsDataFrame.constData() +
                                          (offset * 4) + (i * 8) + 4));
    }
  qDebug();
  return scan_intensities;
}
 
 
void
TimsFrame::cumulateScan(std::size_t scanNum,
                        std::map<quint32, quint32> &accumulate_into) const
{
  qDebug();
 
  if(m_timsDataFrame.size() == 0)
    return;
  checkScanNum(scanNum);
 
 
  std::size_t size = getNbrPeaks(scanNum);
 
  std::size_t offset = getScanOffset(scanNum);
 
  qint32 previous = -1;
  for(std::size_t i = 0; i < size; i++)
    {
      quint32 x =
        (*(quint32 *)((m_timsDataFrame.constData() + (offset * 4) + (i * 8))) +
         previous);
      quint32 y = (*(quint32 *)(m_timsDataFrame.constData() + (offset * 4) +
                                (i * 8) + 4));
 
      previous = x;
 
      auto ret = accumulate_into.insert(std::pair<quint32, quint32>(x, y));
 
      if(ret.second == false)
        {
          // already existed : cumulate
          ret.first->second += y;
        }
    }
  qDebug();
}
 
 
Trace
TimsFrame::cumulateScanToTrace(std::size_t scanNumBegin,
                               std::size_t scanNumEnd) const
{
  qDebug();
  try
    {
      Trace new_trace;
      if(m_timsDataFrame.size() == 0)
        return new_trace;
      std::map<quint32, quint32> raw_spectrum;
      // double local_accumulationTime = 0;
 
      std::size_t imax = scanNumEnd + 1;
      qDebug();
      for(std::size_t i = scanNumBegin; i < imax; i++)
        {
          qDebug() << i;
          cumulateScan(i, raw_spectrum);
          qDebug() << i;
          // local_accumulationTime += m_accumulationTime;
        }
      qDebug();
      pappso::DataPoint data_point_cumul;
 
      for(std::pair<quint32, quint32> pair_tof_intensity : raw_spectrum)
        {
          data_point_cumul.x =
            getMzFromTof(getTofFromIndex(pair_tof_intensity.first));
          // normalization
          data_point_cumul.y =
            pair_tof_intensity.second * ((double)100.0 / m_accumulationTime);
          new_trace.push_back(data_point_cumul);
        }
      new_trace.sortX();
      qDebug();
      return new_trace;
    }
 
  catch(std::exception &error)
    {
      qDebug() << __FILE__ << "@" << __LINE__ << __FUNCTION__ << "()"
               << QString(
                    "Failure in TimsFrame::cumulateScanToTrace %1 to %2 :\n %3")
                    .arg(scanNumBegin, scanNumEnd)
                    .arg(error.what());
    }
}
 
 
pappso::Trace
TimsFrame::getTraceFromCumulatedScans(
  std::map<quint32, quint32> &accumulated_scans) const
{
  qDebug();
  // qDebug();
  // add flanking peaks
  std::vector<quint32> keys;
  transform(begin(accumulated_scans),
            end(accumulated_scans),
            back_inserter(keys),
            [](std::map<quint32, quint32>::value_type const &pair) {
              return pair.first;
            });
  std::sort(keys.begin(), keys.end());
  pappso::DataPoint data_point_cumul;
  data_point_cumul.x = 0;
  quint32 previous   = 0;
  quint32 hole_width = 2;
 
  pappso::Trace local_trace;
 
  bool go_up                = true;
  double previous_intensity = 0;
  /*
  for(quint32 key : keys)
  {
 
              data_point_cumul.x =
                getMzFromTof(getTofFromIndex(key));
              // normalization
              data_point_cumul.y =
                accumulated_scans[key] * ((double)100.0 / m_accumulationTime);
              mass_spectrum_sptr.get()->push_back(data_point_cumul);
    }
*/
  for(quint32 key : keys)
    {
      if((key - previous) > hole_width)
        {
 
          // qDebug() << " key:" << key << " previous:" << previous
          //       << " (key - previous):" << (key - previous);
          if(data_point_cumul.x != 0)
            {
              data_point_cumul.x = data_point_cumul.x / data_point_cumul.y;
 
              data_point_cumul.x =
                getMzFromTof(getTofFromIndex(data_point_cumul.x));
              // normalization
              data_point_cumul.y =
                data_point_cumul.y * ((double)100.0 / m_accumulationTime);
              local_trace.push_back(data_point_cumul);
            }
          // reset
          data_point_cumul.x = 0;
          data_point_cumul.y = 0;
          previous_intensity = 0;
          go_up              = true;
        }
      double intensity = accumulated_scans[key];
 
      if(previous_intensity < intensity)
        {
          if(!go_up)
            {
              // stop
              data_point_cumul.x = data_point_cumul.x / data_point_cumul.y;
 
              data_point_cumul.x =
                getMzFromTof(getTofFromIndex(data_point_cumul.x));
              // normalization
              data_point_cumul.y =
                data_point_cumul.y * ((double)100.0 / m_accumulationTime);
              local_trace.push_back(data_point_cumul);
 
              // reset
              data_point_cumul.x = 0;
              data_point_cumul.y = 0;
              previous_intensity = 0;
              go_up              = true;
            }
        }
      else
        {
          if(go_up)
            {
              go_up = false;
            }
          else
            {
            }
        }
 
 
      data_point_cumul.x += ((double)key * intensity);
      data_point_cumul.y += intensity;
      previous           = key;
      previous_intensity = intensity;
    }
 
 
  if(data_point_cumul.x != 0)
    {
      data_point_cumul.x = data_point_cumul.x / data_point_cumul.y;
 
      data_point_cumul.x = getMzFromTof(getTofFromIndex(data_point_cumul.x));
      // normalization
      data_point_cumul.y =
        data_point_cumul.y * ((double)100.0 / m_accumulationTime);
      local_trace.push_back(data_point_cumul);
    }
  local_trace.sortX();
  qDebug();
  // qDebug();
  return local_trace;
}
 
pappso::MassSpectrumCstSPtr
TimsFrame::getMassSpectrumCstSPtr(std::size_t scanNum) const
{
  qDebug();
  return getMassSpectrumSPtr(scanNum);
}
 
pappso::MassSpectrumSPtr
TimsFrame::getMassSpectrumSPtr(std::size_t scanNum) const
{
 
  qDebug() << " scanNum=" << scanNum;
 
  checkScanNum(scanNum);
 
  qDebug();
 
  pappso::MassSpectrumSPtr mass_spectrum_sptr =
    std::make_shared<pappso::MassSpectrum>();
  // std::vector<DataPoint>
 
  if(m_timsDataFrame.size() == 0)
    return mass_spectrum_sptr;
  qDebug();
 
  std::size_t size = getNbrPeaks(scanNum);
 
  std::size_t offset = getScanOffset(scanNum);
 
  qint32 previous = -1;
  std::vector<quint32> index_list;
  for(std::size_t i = 0; i < size; i++)
    {
      DataPoint data_point(
        (*(quint32 *)((m_timsDataFrame.constData() + (offset * 4) + (i * 8))) +
         previous),
        (*(quint32 *)(m_timsDataFrame.constData() + (offset * 4) + (i * 8) +
                      4)));
 
      // intensity normalization
      data_point.y *= 100.0 / m_accumulationTime;
 
      previous = data_point.x;
 
 
      // mz calibration
      double tof   = (data_point.x * m_digitizerTimebase) + m_digitizerDelay;
      data_point.x = getMzFromTof(tof);
      mass_spectrum_sptr.get()->push_back(data_point);
    }
  qDebug();
  return mass_spectrum_sptr;
}
 
 
void
TimsFrame::extractTimsXicListInRtRange(
  std::vector<TimsXicStructure>::iterator &itXicListbegin,
  std::vector<TimsXicStructure>::iterator &itXicListend,
  XicExtractMethod method) const
{
  std::vector<TimsXicStructure> tmp_xic_list(itXicListbegin, itXicListend);
 
  if(tmp_xic_list.size() == 0)
    return;
  /*
  std::sort(tmp_xic_list.begin(), tmp_xic_list.end(), [](const TimsXicStructure
  &a, const TimsXicStructure &b) { return a.mobilityIndexBegin <
  b.mobilityIndexBegin;
            });
            */
  std::vector<std::size_t> unique_scan_num_list;
  for(auto &&struct_xic : tmp_xic_list)
    {
      for(std::size_t scan = struct_xic.mobilityIndexBegin;
          scan <= struct_xic.mobilityIndexEnd;
          scan++)
        {
          unique_scan_num_list.push_back(scan);
        }
    }
  std::sort(unique_scan_num_list.begin(), unique_scan_num_list.end());
  auto it_scan_num_end =
    std::unique(unique_scan_num_list.begin(), unique_scan_num_list.end());
  auto it_scan_num = unique_scan_num_list.begin();
 
  while(it_scan_num != it_scan_num_end)
    {
      MassSpectrumCstSPtr ms_spectrum = getMassSpectrumCstSPtr(*it_scan_num);
      for(auto &&tmp_xic_struct : tmp_xic_list)
        {
          if(((*it_scan_num) >= tmp_xic_struct.mobilityIndexBegin) &&
             ((*it_scan_num) <= tmp_xic_struct.mobilityIndexEnd))
            {
              if(method == XicExtractMethod::max)
                {
                  tmp_xic_struct.tmpIntensity +=
                    ms_spectrum.get()->maxY(tmp_xic_struct.mzRange.lower(),
                                            tmp_xic_struct.mzRange.upper());
                }
              else
                {
                  // sum
                  tmp_xic_struct.tmpIntensity +=
                    ms_spectrum.get()->sumY(tmp_xic_struct.mzRange.lower(),
                                            tmp_xic_struct.mzRange.upper());
                }
            }
        }
      it_scan_num++;
    }
 
  for(auto &&tmp_xic_struct : tmp_xic_list)
    {
      if(tmp_xic_struct.tmpIntensity != 0)
        {
          tmp_xic_struct.xicSptr.get()->push_back(
            {m_time, tmp_xic_struct.tmpIntensity});
        }
    }
}
 
} // namespace pappso