msrunxicextractordisk.cpp

Go to the documentation of this file.

/**
 * \file pappsomspp/xicextractor/private/msrunxicextractordisk.cpp
 * \date 12/05/2018
 * \author Olivier Langella
 * \brief proteowizard based XIC extractor featuring disk cache
 */
 
/*******************************************************************************
 * Copyright (c) 2018 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 "msrunxicextractordisk.h"
#include <QDebug>
#include "pappsomspp/core/pappsoexception.h"
#include "pappsomspp/core/massspectrum/massspectrum.h"
 
namespace pappso
{
 
MsRunXicExtractorDisk::MsRunXicExtractorDisk(MsRunReaderSPtr &msrun_reader,
                                             const QDir &temporary_dir)
  : pappso::MsRunXicExtractor(msrun_reader)
{
  mpa_temporaryDirectory = nullptr;
  m_temporaryDirectory   = temporary_dir.absolutePath();
}
 
MsRunXicExtractorDisk::MsRunXicExtractorDisk(const MsRunXicExtractorDisk &other)
  : pappso::MsRunXicExtractor(other)
{
 
  m_temporaryDirectory = other.m_temporaryDirectory;
  mpa_temporaryDirectory =
    new QTemporaryDir(QString("%1/msrun_%2_")
                        .arg(m_temporaryDirectory)
                        .arg(msp_msrun_reader.get()->getMsRunId().get()->getXmlId()));
}
 
MsRunXicExtractorDisk::~MsRunXicExtractorDisk()
{
  if(mpa_temporaryDirectory != nullptr)
    {
      delete mpa_temporaryDirectory;
    }
}
 
void
MsRunXicExtractorDisk::prepareExtractor()
{
  qDebug();
  try
    {
      serializeMsRun();
      // msp_msrun_reader = nullptr;
    }
  catch(pappso::PappsoException &errora)
    {
      qDebug();
      throw pappso::PappsoException(QObject::tr("Error reading file (%1) : %2")
                                      .arg(msp_msrun_reader.get()->getMsRunId().get()->toString())
                                      .arg(errora.qwhat()));
    }
  catch(std::exception &error)
    {
      qDebug();
      throw pappso::PappsoException(QObject::tr("Error reading file (%1) using : %2")
                                      .arg(msp_msrun_reader.get()->getMsRunId().get()->toString())
                                      .arg(error.what()));
    }
}
 
 
void
MsRunXicExtractorDisk::protectedExtractXicCoordSPtrList(
  UiMonitorInterface &monitor,
  std::vector<XicCoordSPtr>::iterator it_xic_coord_list_begin,
  std::vector<XicCoordSPtr>::iterator it_xic_coord_list_end)
{
  // sort xic by mz:
  std::sort(it_xic_coord_list_begin, it_xic_coord_list_end, [](XicCoordSPtr &a, XicCoordSPtr &b) {
    return a.get()->mzRange.getMz() < b.get()->mzRange.getMz();
  });
 
  for(auto it = it_xic_coord_list_begin; it != it_xic_coord_list_end; it++)
    {
      extractOneXicCoord(*(it->get()));
      monitor.count();
    }
}
 
void
MsRunXicExtractorDisk::extractOneXicCoord(XicCoord &xic_coord)
{
  std::shared_ptr<Xic> msrunxic_sp = xic_coord.xicSptr;
 
  double rt_begin = xic_coord.rtTarget - m_retentionTimeAroundTarget;
  double rt_end   = xic_coord.rtTarget + m_retentionTimeAroundTarget;
 
 
  std::vector<MsRunSliceSPtr> slice_list;
  slice_list = acquireSlices(xic_coord.mzRange);
 
  if(slice_list.size() == 0)
    {
      throw pappso::PappsoException(QObject::tr("Error getMsRunXicSp slice_list.size() == 0"));
    }
 
  for(std::size_t i = 0; i < m_retentionTimeList.size(); i++)
    {
 
      DataPoint xic_element;
      xic_element.x = m_retentionTimeList[i];
      xic_element.y = 0;
      if((xic_element.x < rt_begin) || (xic_element.x > rt_end))
        continue;
 
      for(auto &&msrun_slice : slice_list)
        {
          const MassSpectrum &spectrum = msrun_slice.get()->getSpectrum(i);
          for(auto &&peak : spectrum)
            {
              if(xic_coord.mzRange.contains(peak.x))
                {
                  if(m_xicExtractMethod == Enums::XicExtractMethod::sum)
                    {
                      xic_element.y += peak.y;
                    }
                  else
                    {
                      if(xic_element.y < peak.y)
                        {
                          xic_element.y = peak.y;
                        }
                    }
                }
            }
        }
      msrunxic_sp.get()->push_back(xic_element);
    }
}
 
void
MsRunXicExtractorDisk::serializeMsRun()
{
  qDebug();
  m_minMz = 5000;
  m_maxMz = 0;
 
  unsigned int slice_number;
  std::map<unsigned int, MassSpectrum> spectrum_map;
 
  /*
    const pwiz::msdata::SpectrumList *p_spectrum_list =
      p_msdatafile->run.spectrumListPtr.get();
 
    std::size_t spectrum_list_size = p_spectrum_list->size();
    pwiz::msdata::SpectrumPtr pwiz_spectrum;
    */
 
  m_rtSize = m_msrun_points.size();
 
 
  MassSpectrumCstSPtr spectrum;
  for(auto &&msrun_point : m_msrun_points)
    {
 
      spectrum_map.clear();
 
      m_retentionTimeList.push_back(msrun_point.rt);
 
      spectrum = msp_msrun_reader.get()->massSpectrumCstSPtr(msrun_point.spectrum_index);
 
      const MassSpectrum *p_spectrum = spectrum.get();
      if(p_spectrum->size() > 0)
        {
          if(p_spectrum->begin()->x < m_minMz)
            {
              m_minMz = p_spectrum->begin()->x;
            }
          // iterate through the m/z-intensity pairs
 
          if(p_spectrum->back().x > m_maxMz)
            {
              m_maxMz = p_spectrum->back().x;
            }
 
          for(auto &peak : *p_spectrum)
            {
 
              slice_number = peak.x;
 
              std::pair<std::map<unsigned int, MassSpectrum>::iterator, bool> ret =
                spectrum_map.insert(
                  std::pair<unsigned int, MassSpectrum>(slice_number, MassSpectrum()));
 
              ret.first->second.push_back(peak);
              // auto ret = spectrum_map.insert(std::pair<unsigned int,
              // MassSpectrum>(slice_number,MassSpectrum()));
              // ret.first->second.push_back(peak);
            }
 
          // slices are ready for this retention time
          storeSlices(spectrum_map, m_retentionTimeList.size() - 1);
        }
    }
 
  endPwizRead();
  qDebug();
}
 
 
void
MsRunXicExtractorDisk::storeSlices(std::map<unsigned int, MassSpectrum> &spectrum_map,
                                   std::size_t ipos)
{
  // qDebug() << __FILE__ << " " << __FUNCTION__ << " " << __LINE__;
 
  for(auto &&spectrum_pair : spectrum_map)
    {
      appendSliceOnDisk(spectrum_pair.first, spectrum_pair.second, ipos);
    }
 
  // qDebug() << __FILE__ << " " << __FUNCTION__ << " " << __LINE__;
}
 
void
MsRunXicExtractorDisk::appendSliceOnDisk(unsigned int slice_number,
                                         MassSpectrum &spectrum,
                                         std::size_t ipos)
{
  // qDebug() << __FILE__ << " " << __FUNCTION__ << " " << __LINE__;
  QFile slice_file(QString("%1/%2").arg(mpa_temporaryDirectory->path()).arg(slice_number));
  bool new_file = false;
  if(!slice_file.exists())
    {
      new_file = true;
    }
  if(!slice_file.open(QIODevice::WriteOnly | QIODevice::Append))
    {
      throw pappso::PappsoException(
        QObject::tr("unable to open file %1").arg(slice_file.fileName()));
    }
  QDataStream stream(&slice_file);
 
  if(new_file)
    {
      stream << (quint32)slice_number;
      stream << (quint32)m_rtSize;
    }
 
  stream << (quint32)ipos;
  stream << spectrum;
 
  slice_file.flush();
  slice_file.close();
  // qDebug() << __FILE__ << " " << __FUNCTION__ << " " << __LINE__;
}
 
MsRunSliceSPtr
MsRunXicExtractorDisk::unserializeSlice(unsigned int slice_number)
{
  qDebug();
  try
    {
      std::shared_ptr<MsRunSlice> msrun_slice_sp = std::make_shared<MsRunSlice>(MsRunSlice());
 
      QFile slice_file(QString("%1/%2").arg(mpa_temporaryDirectory->path()).arg(slice_number));
      if(!slice_file.exists())
        {
          msrun_slice_sp.get()->setSize(m_rtSize);
          msrun_slice_sp.get()->setSliceNumber(slice_number);
          return msrun_slice_sp;
        }
      if(!slice_file.open(QIODevice::ReadOnly))
        {
          throw pappso::PappsoException(
            QObject::tr("unable to open file %1 in readonly").arg(slice_file.fileName()));
        }
      QDataStream stream(&slice_file);
 
      stream >> *(msrun_slice_sp.get());
 
      slice_file.close();
 
      return msrun_slice_sp;
    }
  catch(pappso::PappsoException &error)
    {
      throw pappso::PappsoException(
        QObject::tr("error unserializing slice %1:\n%2").arg(slice_number).arg(error.qwhat()));
    }
  qDebug();
}
 
std::vector<MsRunSliceSPtr>
MsRunXicExtractorDisk::acquireSlices(const MzRange &mz_range)
{
  QMutexLocker lock(&m_mutex);
  std::vector<MsRunSliceSPtr> slice_list;
  for(unsigned int i = mz_range.lower(); i <= mz_range.upper(); i++)
    {
      auto it = std::find_if(
        m_msRunSliceListCache.begin(),
        m_msRunSliceListCache.end(),
        [i](const MsRunSliceSPtr &slice_sp) { return slice_sp.get()->getSliceNumber() == i; });
      if(it != m_msRunSliceListCache.end())
        {
          slice_list.push_back(*it);
          m_msRunSliceListCache.push_back(*it);
        }
      else
        {
          MsRunSliceSPtr slice_sp = unserializeSlice(i);
          slice_list.push_back(slice_sp);
          m_msRunSliceListCache.push_back(slice_sp);
        }
    }
 
  if(m_msRunSliceListCache.size() > 20)
    {
      m_msRunSliceListCache.pop_front();
    }
  return slice_list;
}
 
 
void
MsRunXicExtractorDisk::endPwizRead()
{
  msp_msrun_reader.get()->releaseDevice();
}
} // namespace pappso