OpenSwathBase.h

Go to the documentation of this file.

// --------------------------------------------------------------------------
//                   OpenMS -- Open-Source Mass Spectrometry
// --------------------------------------------------------------------------
// Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
// ETH Zurich, and Freie Universitaet Berlin 2002-2017.
//
// This software is released under a three-clause BSD license:
//  * Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//  * Redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//  * Neither the name of any author or any participating institution
//    may be used to endorse or promote products derived from this software
//    without specific prior written permission.
// For a full list of authors, refer to the file AUTHORS.
// --------------------------------------------------------------------------
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING
// INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
// OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
// WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
// ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// --------------------------------------------------------------------------
// $Maintainer: Hannes Roest$
// $Authors: Hannes Roest$
// --------------------------------------------------------------------------

#pragma once

// Consumers
#include <OpenMS/FORMAT/DATAACCESS/MSDataWritingConsumer.h>
#include <OpenMS/FORMAT/DATAACCESS/MSDataSqlConsumer.h>

// Files
#include <OpenMS/FORMAT/FileHandler.h>
#include <OpenMS/FORMAT/FileTypes.h>
#include <OpenMS/FORMAT/TraMLFile.h>
#include <OpenMS/FORMAT/MzMLFile.h>
#include <OpenMS/FORMAT/FeatureXMLFile.h>
#include <OpenMS/FORMAT/TransformationXMLFile.h>
#include <OpenMS/FORMAT/SwathFile.h>
#include <OpenMS/ANALYSIS/OPENSWATH/SwathWindowLoader.h>
#include <OpenMS/ANALYSIS/OPENSWATH/TransitionTSVFile.h>
#include <OpenMS/ANALYSIS/OPENSWATH/TransitionPQPFile.h>
#include <OpenMS/ANALYSIS/OPENSWATH/OpenSwathTSVWriter.h>
#include <OpenMS/ANALYSIS/OPENSWATH/OpenSwathOSWWriter.h>

// Kernel and implementations
#include <OpenMS/KERNEL/MSExperiment.h>
#include <OpenMS/ANALYSIS/OPENSWATH/DATAACCESS/SpectrumAccessOpenMS.h>
#include <OpenMS/ANALYSIS/OPENSWATH/DATAACCESS/SpectrumAccessTransforming.h>
#include <OpenMS/ANALYSIS/OPENSWATH/DATAACCESS/SpectrumAccessOpenMSInMemory.h>
#include <OpenMS/OPENSWATHALGO/DATAACCESS/SwathMap.h>

// Helpers
#include <OpenMS/ANALYSIS/OPENSWATH/OpenSwathHelper.h>
#include <OpenMS/ANALYSIS/OPENSWATH/DATAACCESS/DataAccessHelper.h>
#include <OpenMS/ANALYSIS/OPENSWATH/DATAACCESS/SimpleOpenMSSpectraAccessFactory.h>

// Algorithms
#include <OpenMS/ANALYSIS/OPENSWATH/MRMRTNormalizer.h>
#include <OpenMS/ANALYSIS/OPENSWATH/ChromatogramExtractor.h>
#include <OpenMS/ANALYSIS/OPENSWATH/MRMFeatureFinderScoring.h>
#include <OpenMS/ANALYSIS/OPENSWATH/MRMTransitionGroupPicker.h>
#include <OpenMS/ANALYSIS/OPENSWATH/SwathMapMassCorrection.h>

#include <OpenMS/ANALYSIS/OPENSWATH/OpenSwathWorkflow.h>

#include <cassert>
#include <limits>


#include <OpenMS/APPLICATIONS/TOPPBase.h>



namespace OpenMS
{

class TOPPOpenSwathBase :
  public TOPPBase
{

public:

  TOPPOpenSwathBase(String name, String description, bool official = true) :
    TOPPBase(name, description, official)
  {
  }

private:

  void loadSwathFiles_(const StringList& file_list,
                       const bool split_file,
                       const String& tmp,
                       const String& readoptions,
                       boost::shared_ptr<ExperimentalSettings > & exp_meta,
                       std::vector< OpenSwath::SwathMap > & swath_maps,
                       Interfaces::IMSDataConsumer* plugin_consumer)
  {
    SwathFile swath_file;
    swath_file.setLogType(log_type_);

    if (split_file || file_list.size() > 1)
    {
      // TODO cannot use data reduction here any more ...
      swath_maps = swath_file.loadSplit(file_list, tmp, exp_meta, readoptions);
    }
    else
    {
      FileTypes::Type in_file_type = FileHandler::getTypeByFileName(file_list[0]);
      if (in_file_type == FileTypes::MZML)
      {
        swath_maps = swath_file.loadMzML(file_list[0], tmp, exp_meta, readoptions, plugin_consumer);
      }
      else if (in_file_type == FileTypes::MZXML)
      {
        swath_maps = swath_file.loadMzXML(file_list[0], tmp, exp_meta, readoptions);
      }
      else if (in_file_type == FileTypes::SQMASS)
      {
        swath_maps = swath_file.loadSqMass(file_list[0], exp_meta);
      }
      else
      {
        throw Exception::IllegalArgument(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION,
            "Input file needs to have ending mzML or mzXML");
      }
    }
  }

protected:

  bool loadSwathFiles(const StringList& file_list,
                      boost::shared_ptr<ExperimentalSettings >& exp_meta,
                      std::vector< OpenSwath::SwathMap >& swath_maps,
                      const bool split_file,
                      const String& tmp,
                      const String& readoptions,
                      const String& swath_windows_file,
                      const double min_upper_edge_dist,
                      const bool force,
                      const bool sort_swath_maps,
                      const bool sonar,
                      Interfaces::IMSDataConsumer* plugin_consumer = nullptr)
  {
    // (i) Load files
    loadSwathFiles_(file_list, split_file, tmp, readoptions, exp_meta, swath_maps, plugin_consumer);

    // (ii) Allow the user to specify the SWATH windows
    if (!swath_windows_file.empty())
    {
      SwathWindowLoader::annotateSwathMapsFromFile(swath_windows_file, swath_maps, sort_swath_maps, force);
    }

    for (Size i = 0; i < swath_maps.size(); i++)
    {
      LOG_DEBUG << "Found swath map " << i
        << " with lower " << swath_maps[i].lower
        << " and upper " << swath_maps[i].upper
        << " and " << swath_maps[i].sptr->getNrSpectra()
        << " spectra." << std::endl;
    }

    // (iii) Sanity check: there should be no overlap between the windows:
    std::vector<std::pair<double, double>> sw_windows;
    for (Size i = 0; i < swath_maps.size(); i++)
    {
      if (!swath_maps[i].ms1)
      {
        sw_windows.push_back(std::make_pair(swath_maps[i].lower, swath_maps[i].upper));
      }
    }
    // sort by lower bound (first entry in pair)
    std::sort(sw_windows.begin(), sw_windows.end());

    for (Size i = 1; i < sw_windows.size(); i++)
    {
      double lower_map_end = sw_windows[i-1].second - min_upper_edge_dist;
      double upper_map_start = sw_windows[i].first;
      LOG_DEBUG << "Extraction will go up to " << lower_map_end << " and continue at " << upper_map_start << std::endl;

      if (upper_map_start - lower_map_end > 0.01)
      {
        LOG_WARN << "Extraction will have a gap between " << lower_map_end << " and " << upper_map_start << std::endl;
        if (!force)
        {
          LOG_ERROR << "Extraction windows have a gap. Will abort (override with -force)" << std::endl;
          return false;
        }
      }

      if (sonar) {continue;} // skip next step as expect them to overlap ...

      if (lower_map_end - upper_map_start > 0.01)
      {
        LOG_WARN << "Extraction will overlap between " << lower_map_end << " and " << upper_map_start << "!\n"
                 << "This will lead to multiple extraction of the transitions in the overlapping region "
                 << "which will lead to duplicated output. It is very unlikely that you want this." << "\n"
                 << "Please fix this by providing an appropriate extraction file with -swath_windows_file" << std::endl;
        if (!force)
        {
          LOG_ERROR << "Extraction windows overlap. Will abort (override with -force)" << std::endl;
          return false;
        }
      }
    }
    return true;
  }

  void prepareChromOutput(Interfaces::IMSDataConsumer ** chromatogramConsumer, 
                          const boost::shared_ptr<ExperimentalSettings>& exp_meta,
                          const OpenSwath::LightTargetedExperiment& transition_exp,
                          const String& out_chrom)
  {
    if (!out_chrom.empty())
    {
      if (out_chrom.hasSuffix(".sqMass"))
      {
        bool full_meta = false; // can lead to very large files in memory
        bool lossy_compression = true;
        *chromatogramConsumer = new MSDataSqlConsumer(out_chrom, 500, full_meta, lossy_compression);
      }
      else
      {
        PlainMSDataWritingConsumer * chromConsumer = new PlainMSDataWritingConsumer(out_chrom);
        int expected_chromatograms = transition_exp.transitions.size();
        chromConsumer->setExpectedSize(0, expected_chromatograms);
        chromConsumer->setExperimentalSettings(*exp_meta);
        chromConsumer->getOptions().setWriteIndex(true);  // ensure that we write the index
        chromConsumer->addDataProcessing(getProcessingInfo_(DataProcessing::SMOOTHING));

        // prepare data structures for lossy compression
        MSNumpressCoder::NumpressConfig npconfig_mz;
        MSNumpressCoder::NumpressConfig npconfig_int;
        npconfig_mz.estimate_fixed_point = true; // critical
        npconfig_int.estimate_fixed_point = true; // critical
        npconfig_mz.numpressErrorTolerance = -1.0; // skip check, faster
        npconfig_int.numpressErrorTolerance = -1.0; // skip check, faster
        npconfig_mz.setCompression("linear");
        npconfig_int.setCompression("slof");
        npconfig_mz.linear_fp_mass_acc = 0.05; // set the desired RT accuracy in seconds

        chromConsumer->getOptions().setNumpressConfigurationMassTime(npconfig_mz);
        chromConsumer->getOptions().setNumpressConfigurationIntensity(npconfig_int);
        chromConsumer->getOptions().setCompression(true);

        *chromatogramConsumer = chromConsumer;
      }
    }
    else
    {
      *chromatogramConsumer = new NoopMSDataWritingConsumer("");
    }
  }

  OpenSwath::LightTargetedExperiment loadTransitionList(const FileTypes::Type& tr_type,
                                                        const String& tr_file,
                                                        const Param& tsv_reader_param)
  {
    OpenSwath::LightTargetedExperiment transition_exp;
    ProgressLogger progresslogger;
    progresslogger.setLogType(log_type_);
    if (tr_type == FileTypes::TRAML)
    {
      progresslogger.startProgress(0, 1, "Load TraML file");
      TargetedExperiment targeted_exp;
      TraMLFile().load(tr_file, targeted_exp);
      OpenSwathDataAccessHelper::convertTargetedExp(targeted_exp, transition_exp);
      progresslogger.endProgress();
    }
    else if (tr_type == FileTypes::PQP)
    {
      progresslogger.startProgress(0, 1, "Load PQP file");
      TransitionPQPFile().convertPQPToTargetedExperiment(tr_file.c_str(), transition_exp);
      progresslogger.endProgress();
    }
    else if (tr_type == FileTypes::TSV)
    {
      progresslogger.startProgress(0, 1, "Load TSV file");
      TransitionTSVFile tsv_reader;
      tsv_reader.setParameters(tsv_reader_param);
      tsv_reader.convertTSVToTargetedExperiment(tr_file.c_str(), tr_type, transition_exp);
      progresslogger.endProgress();
    }
    else
    {
      LOG_ERROR << "Provide valid TraML, TSV or PQP transition file." << std::endl;
      throw Exception::IllegalArgument(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "Need to provide valid input file.");
    }
    return transition_exp;
  }

  TransformationDescription performCalibration(String trafo_in,
        String irt_tr_file,
        std::vector< OpenSwath::SwathMap > & swath_maps,
        double min_rsq,
        double min_coverage,
        const Param& feature_finder_param,
        const ChromExtractParams& cp_irt,
        const Param& irt_detection_param,
        const String & mz_correction_function,
        Size debug_level,
        bool sonar,
        bool load_into_memory,
        const String& irt_trafo_out,
        const String& irt_mzml_out)
  {
    TransformationDescription trafo_rtnorm;

    if (!trafo_in.empty())
    {
      // get read RT normalization file
      TransformationXMLFile trafoxml;
      trafoxml.load(trafo_in, trafo_rtnorm, false);
      Param model_params = getParam_().copy("model:", true);
      model_params.setValue("symmetric_regression", "false");
      model_params.setValue("span", irt_detection_param.getValue("lowess:span"));
      model_params.setValue("num_nodes", irt_detection_param.getValue("b_spline:num_nodes"));
      String model_type = irt_detection_param.getValue("alignmentMethod");
      trafo_rtnorm.fitModel(model_type, model_params);
    }
    else if (!irt_tr_file.empty())
    {
      // Loading iRT file
      std::cout << "Will load iRT transitions and try to find iRT peptides" << std::endl;
      TraMLFile traml;
      FileTypes::Type tr_type = FileHandler::getType(irt_tr_file);
      Param tsv_reader_param = TransitionTSVFile().getDefaults();
      OpenSwath::LightTargetedExperiment irt_transitions = loadTransitionList(tr_type, irt_tr_file, tsv_reader_param);

      // perform extraction
      OpenSwathCalibrationWorkflow wf;
      wf.setLogType(log_type_);
      trafo_rtnorm = wf.performRTNormalization(irt_transitions, swath_maps, min_rsq, min_coverage,
      feature_finder_param, cp_irt, irt_detection_param, mz_correction_function, irt_mzml_out,
      debug_level, sonar, load_into_memory);

      if (!irt_trafo_out.empty())
      {
        TransformationXMLFile().store(irt_trafo_out, trafo_rtnorm);
      }
    }
    return trafo_rtnorm;
  }


};

}