mzxmloutput.cpp

Go to the documentation of this file.

/**
 * \file pappsomspp/msrun/output/mzxmloutput.cpp
 * \date 23/11/2019
 * \author Olivier Langella
 * \brief write msrun peaks into mzxml output stream
 */
 
/*******************************************************************************
 * 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 "mzxmloutput.h"
#include <QDebug>
#include <QStringList>
#include <algorithm>
#include <cstdio>
#include "pappsomspp/config.h"
 
using namespace pappso;
 
 
template <class T>
T
change_endian(T in)
{
  char *const p = reinterpret_cast<char *>(&in);
  for(size_t i = 0; i < sizeof(T) / 2; ++i)
    std::swap(p[i], p[sizeof(T) - i - 1]);
  return in;
}
 
 
MzxmlOutput::Translater::Translater(MzxmlOutput *p_mzxml_output)
{
  mp_output = p_mzxml_output;
}
MzxmlOutput::Translater::~Translater()
{
}
void
MzxmlOutput::Translater::setQualifiedMassSpectrum(const QualifiedMassSpectrum &spectrum)
{
  qDebug();
  mp_output->m_monitor.count();
  mp_output->writeQualifiedMassSpectrum(spectrum);
  qDebug();
}
bool
MzxmlOutput::Translater::needPeakList() const
{
  return true;
}
 
 
MzxmlOutput::MzxmlOutput(UiMonitorInterface &monitor, QIODevice *p_output_device)
  : m_monitor(monitor)
{
 
  mpa_outputStream = new QXmlStreamWriter(p_output_device);
  mpa_outputStream->setAutoFormatting(true);
 
  mpa_outputStream->writeStartDocument("1.0");
}
 
MzxmlOutput::~MzxmlOutput()
{
  close();
  delete mpa_outputStream;
}
 
void
MzxmlOutput::setReadAhead(bool isReadAhead)
{
  m_isReadAhead = isReadAhead;
}
void
MzxmlOutput::write(MsRunReader *p_msrunreader)
{
  qDebug();
  m_monitor.setTotalSteps(p_msrunreader->spectrumListSize());
  writeHeader(p_msrunreader);
 
  Translater translater(this);
 
  translater.setReadAhead(m_isReadAhead);
 
  translater.setNeedMsLevelPeakList(1, !m_ms1IsMasked);
  // translater.setNeedMsLevelPeakList(1, false);
  // translater.setNeedMsLevelPeakList(2, false);
 
 
  MsRunReadConfig config;
  std::vector<size_t> ms_levels;
  for(std::size_t i = 1; i < 9; i++)
    {
      ms_levels.push_back(i);
    }
  config.setMsLevels(ms_levels);
  config.setNeedPeakList(true);
 
  p_msrunreader->readSpectrumCollection2(config, translater);
 
  m_monitor.setTotalSteps(0);
  qDebug();
}
 
void
pappso::MzxmlOutput::write(const MsRunReadConfig &read_config, pappso::MsRunReader *p_msrunreader)
{
  qDebug();
  m_monitor.setTotalSteps(p_msrunreader->spectrumListSize());
  writeHeader(p_msrunreader);
 
  Translater translater(this);
 
  translater.setReadAhead(m_isReadAhead);
 
  translater.setNeedMsLevelPeakList(1, !m_ms1IsMasked);
  // translater.setNeedMsLevelPeakList(1, false);
  // translater.setNeedMsLevelPeakList(2, false);
  p_msrunreader->readSpectrumCollection2(read_config, translater);
 
  m_monitor.setTotalSteps(0);
}
 
void
MzxmlOutput::writeHeader(MsRunReader *p_msrunreader)
{
 
  mpa_outputStream->writeNamespace("http://www.w3.org/2001/XMLSchema-instance", "xsi");
  // xmlns="http://sashimi.sourceforge.net/schema_revision/mzXML_2.0"
  // xsi:schemaLocation="http://sashimi.sourceforge.net/schema_revision/mzXML_2.0
  // http://sashimi.sourceforge.net/schema_revision/mzXML_2.0/mzXML_idx_2.0.xsd"
  /*
114                 writer.setPrefix("xsi", xmlnsxsi);
115                 writer.setDefaultNamespace(namespaceURI);
mpa_outputStream->writeStartElement("mzXML");
117                 writer.writeNamespace("xsi", xmlnsxsi);
118                 writer.writeDefaultNamespace(namespaceURI);
119
120                 writer.writeAttribute(xmlnsxsi, "schemaLocation",
xsischemaLocation); 121 */
  mpa_outputStream->writeStartElement("mzXML");
  mpa_outputStream->writeAttribute("xmlns",
                                   "http://sashimi.sourceforge.net/schema_revision/mzXML_3.2");
  mpa_outputStream->writeAttribute("xsi:schemaLocation",
                                   "http://sashimi.sourceforge.net/schema_revision/mzXML_3.2 "
                                   "http://sashimi.sourceforge.net/schema_revision/mzXML_3.2/"
                                   "mzXML_idx_3.2.xsd");
 
  mpa_outputStream->writeStartElement("msRun");
  mpa_outputStream->writeAttribute("scanCount",
                                   QString("%1").arg(p_msrunreader->spectrumListSize()));
  //<msRun scanCount="16576" startTime="PT0.292553S" endTime="PT3000.34S">
  // writer.writeAttribute("scanCount",
  // ms_run.getSpectrumCount(this.controller).toString());
 
  /*
   * # < parentFile fileName = #
   * "file://SEQUEST1/raw/vidal/20060411_VIDAL_JEAN_1_PEPCR1_42140.RAW" #
   * fileType = "RAWData" fileSha1 = #
   * "23c1620d4ad3f4f0103b0141b7caec1e8b7eebf5" / >
   */
  mpa_outputStream->writeStartElement("parentFile");
  mpa_outputStream->writeAttribute("fileName", p_msrunreader->getMsRunId()->getFileName());
  mpa_outputStream->writeAttribute("fileType", "RAWData");
  mpa_outputStream->writeEndElement();
  /*
144
145                 MsInstrumentList instrument_list =
ms_run.getMsInstruments(controller); 146                 for (MsInstrument
instrument : instrument_list) { 147 this.write(instrument); 148 }
*/
 
  mpa_outputStream->writeStartElement("msInstrument");
  mpa_outputStream->writeAttribute("msInstrumentID", "1");
  //<msManufacturer category="msManufacturer" value="Thermo Scientific"/>
  mpa_outputStream->writeStartElement("msManufacturer");
  mpa_outputStream->writeAttribute("category", "msManufacturer");
  mpa_outputStream->writeAttribute("value", "unknown");
  mpa_outputStream->writeEndElement();
  //<msModel category="msModel" value="Q Exactive"/>
  // <msIonisation category="msIonisation" value="nanoelectrospray"/>
  // <msMassAnalyzer category="msMassAnalyzer" value="quadrupole"/>
  // <msDetector category="msDetector" value="inductive detector"/>
  // <software type="acquisition" name="Xcalibur"
  // version="2.1-152001/2.1.0.1520"/>
  mpa_outputStream->writeEndElement();
  /*
149
150                 // #< dataProcessing centroided ="1" >
151                 // my $ref_data_processings =
$ms_run_description->dataProcessing(); 152                 MsDataProcessingList
dataProcList = ms_run.getMsDataProcessings(controller); 153                 for
(MsDataProcessing msDataProc : dataProcList) { 154 this.write(msDataProc); 155 }
*/
  mpa_outputStream->writeStartElement("dataProcessing");
  //<dataProcessing centroided="1">
  mpa_outputStream->writeAttribute("centroided", "1");
  //  <software type="conversion" name="ProteoWizard" version="3.0.3706"/>
  mpa_outputStream->writeStartElement("software");
  mpa_outputStream->writeAttribute("type", "conversion");
  mpa_outputStream->writeAttribute("name", PAPPSOMSPP_NAME);
  mpa_outputStream->writeAttribute("version", PAPPSOMSPP_VERSION);
  mpa_outputStream->writeEndElement();
  //<processingOperation name="Conversion to mzML"/>
  mpa_outputStream->writeStartElement("processingOperation");
  mpa_outputStream->writeAttribute("name", "Conversion to mzXML");
  //<software type="processing" name="ProteoWizard" version="3.0.3706"/>
  mpa_outputStream->writeStartElement("software");
  mpa_outputStream->writeAttribute("type", "processing");
  mpa_outputStream->writeAttribute("name", PAPPSOMSPP_NAME);
  mpa_outputStream->writeAttribute("version", PAPPSOMSPP_VERSION);
  mpa_outputStream->writeEndElement();
  //<comment>Thermo/Xcalibur peak picking</comment>
  mpa_outputStream->writeStartElement("comment");
  mpa_outputStream->writeCharacters("pappso::MzxmlOutput");
  mpa_outputStream->writeEndElement();
  //</dataProcessing>
  mpa_outputStream->writeEndElement();
  mpa_outputStream->writeEndElement();
  // Peaks
}
 
 
void
MzxmlOutput::close()
{
  mpa_outputStream->writeEndDocument();
}
 
 
std::size_t
MzxmlOutput::getScanNumberFromNativeId(const QString &native_id) const
{
  QStringList native_id_list = native_id.split("=");
  if(native_id_list.size() < 2)
    {
    }
  else
    {
      return native_id_list.back().toULong();
    }
  return std::numeric_limits<std::size_t>::max();
}
 
std::size_t
MzxmlOutput::getScanNumber(const QualifiedMassSpectrum &spectrum) const
{
  std::size_t scan_number = getScanNumberFromNativeId(spectrum.getMassSpectrumId().getNativeId());
  if(scan_number == std::numeric_limits<std::size_t>::max())
    {
      scan_number = spectrum.getMassSpectrumId().getSpectrumIndex() + 1;
    }
  return scan_number;
}
 
std::size_t
MzxmlOutput::getPrecursorScanNumber(const QualifiedMassSpectrum &spectrum) const
{
 
  std::size_t scan_number = getScanNumberFromNativeId(spectrum.getPrecursorNativeId());
  if(scan_number == std::numeric_limits<std::size_t>::max())
    {
      scan_number = spectrum.getPrecursorSpectrumIndex() + 1;
    }
  return scan_number;
}
 
void
MzxmlOutput::writeQualifiedMassSpectrum(const pappso::QualifiedMassSpectrum &spectrum)
{
  qDebug();
  mpa_outputStream->writeStartElement("scan");
  /*
  <scan num="1"
          scanType="Full"
          centroided="1"
          msLevel="1"
          peaksCount="1552"
          polarity="+"
          retentionTime="PT0.292553S"
          lowMz="400.153411865234"
          highMz="1013.123352050781"
          basePeakMz="445.12003"
          basePeakIntensity="2.0422125e06"
          totIonCurrent="1.737798e07">*/
  mpa_outputStream->writeAttribute("num", QString("%1").arg(getScanNumber(spectrum)));
  mpa_outputStream->writeAttribute("centroided", QString("1"));
  mpa_outputStream->writeAttribute("msLevel", QString("%1").arg(spectrum.getMsLevel()));
  if(spectrum.getMassSpectrumCstSPtr().get() == nullptr)
    {
      mpa_outputStream->writeAttribute("peaksCount", "0");
    }
  else
    {
      mpa_outputStream->writeAttribute("peaksCount", QString("%1").arg(spectrum.size()));
 
      if(spectrum.size() > 0)
        {
          mpa_outputStream->writeAttribute(
            "lowMz", QString::number(spectrum.getMassSpectrumCstSPtr().get()->front().x, 'f', 12));
 
          mpa_outputStream->writeAttribute(
            "highMz", QString::number(spectrum.getMassSpectrumCstSPtr().get()->back().x, 'f', 12));
          // mpa_outputStream->writeAttribute("highMz",
          // QString::number(spectrum.getMassSpectrumCstSPtr().get()->back().x,
          // 'f', 10)); basePeakMz="245.1271988"
          //   basePeakIntensity="5810.7739"
          //   totIonCurrent="57803.815999999999">
        }
    }
  mpa_outputStream->writeAttribute("polarity", "+");
  mpa_outputStream->writeAttribute(
    "retentionTime", QString("PT%1S").arg(QString::number(spectrum.getRtInSeconds(), 'f', 2)));
 
  if(spectrum.getMsLevel() > 1)
    {
 
      //<precursorMz precursorScanNum="16574"
      // precursorIntensity="58403.04296875" precursorCharge="2"
      ////activationMethod="HCD">994.690619901808</precursorMz>
      mpa_outputStream->writeStartElement("precursorMz");
      mpa_outputStream->writeAttribute("precursorScanNum",
                                       QString("%1").arg(getPrecursorScanNumber(spectrum)));
      bool ok;
      double precursor_intensity = spectrum.getPrecursorIntensity(&ok);
      if(ok)
        {
          mpa_outputStream->writeAttribute("precursorIntensity",
                                           QString::number(precursor_intensity, 'f', 4));
        }
      uint charge = spectrum.getPrecursorCharge(&ok);
      if(ok)
        {
          mpa_outputStream->writeAttribute("precursorCharge", QString("%1").arg(charge));
        }
      double precursor_mz = spectrum.getPrecursorMz(&ok);
      if(ok)
        {
          mpa_outputStream->writeCharacters(QString::number(precursor_mz, 'f', 12));
        }
      mpa_outputStream->writeEndElement();
    }
 
  /*<peaks compressionType="none"
           compressedLen="0"
           precision="64"
           byteOrder="network"
           contentType="m/z-int"></peaks>*/
 
  mpa_outputStream->writeStartElement("peaks");
  mpa_outputStream->writeAttribute("compressionType", "none");
  mpa_outputStream->writeAttribute("compressedLen", "0");
  mpa_outputStream->writeAttribute("precision", "64");
  mpa_outputStream->writeAttribute("byteOrder", "network");
  mpa_outputStream->writeAttribute("contentType", "m/z-int");
 
  if((spectrum.getMassSpectrumCstSPtr().get() != nullptr) && (!spectrum.isEmptyMassSpectrum()))
    {
      QByteArray byte_array;
      if(QSysInfo::ByteOrder == QSysInfo::LittleEndian)
        {
          for(const DataPoint &peak : *(spectrum.getMassSpectrumCstSPtr().get()))
            {
              double swap = change_endian(peak.x);
              byte_array.append((char *)&swap, 8);
              swap = change_endian(peak.y);
              byte_array.append((char *)&swap, 8);
            }
        }
      else
        {
          for(const DataPoint &peak : *(spectrum.getMassSpectrumCstSPtr().get()))
            {
              byte_array.append((char *)&peak.x, 8);
              byte_array.append((char *)&peak.y, 8);
            }
        }
      mpa_outputStream->writeCharacters(byte_array.toBase64());
    }
  mpa_outputStream->writeEndElement();
 
  // scan
  mpa_outputStream->writeEndElement();
  qDebug();
}
 
void
MzxmlOutput::maskMs1(bool mask_ms1)
{
  m_ms1IsMasked = mask_ms1;
}