FeatureFinderIdentificationAlgorithm.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-2022.
//
// 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: Timo Sachsenberg $
// $Authors: Hendrik Weisser $
// --------------------------------------------------------------------------

#ifndef OPENMS_TRANSFORMATIONS_FEATUREFINDER_FEATUREFINDERIDENTIFICATIONALGORITHM_H
#define OPENMS_TRANSFORMATIONS_FEATUREFINDER_FEATUREFINDERIDENTIFICATIONALGORITHM_H

#include <OpenMS/ANALYSIS/TARGETED/TargetedExperiment.h>
#include <OpenMS/ANALYSIS/MAPMATCHING/TransformationDescription.h>
#include <OpenMS/ANALYSIS/OPENSWATH/MRMFeatureFinderScoring.h>
#include <OpenMS/CONCEPT/ProgressLogger.h>
#include <OpenMS/DATASTRUCTURES/DefaultParamHandler.h>
#include <OpenMS/KERNEL/MSExperiment.h>
#include <OpenMS/TRANSFORMATIONS/FEATUREFINDER/FeatureFinderAlgorithmPickedHelperStructs.h>

#include <vector>
#include <fstream>
#include <map>

namespace OpenMS
{
  class IsotopeDistribution;



class OPENMS_DLLAPI FeatureFinderIdentificationAlgorithm :
  public DefaultParamHandler
{
public:
  FeatureFinderIdentificationAlgorithm(); 

  void run(
    std::vector<PeptideIdentification> peptides,
    const std::vector<ProteinIdentification>& proteins,
    std::vector<PeptideIdentification> peptides_ext,
    std::vector<ProteinIdentification> proteins_ext,
    FeatureMap& features,
    const FeatureMap& seeds = FeatureMap(),
    const String& spectra_file = ""
    );

  void runOnCandidates(FeatureMap& features);

  PeakMap& getMSData();
  const PeakMap& getMSData() const;

  void setMSData(const PeakMap& ms_data); // for pyOpenMS
  void setMSData(PeakMap&& ms_data); // moves peak data and saves the copy. Note that getMSData() will give back a processed/modified version.

  PeakMap& getChromatograms();
  const PeakMap& getChromatograms() const;

  ProgressLogger& getProgressLogger();
  const ProgressLogger& getProgressLogger() const;

  TargetedExperiment& getLibrary();
  const TargetedExperiment& getLibrary() const;

protected:
  typedef FeatureFinderAlgorithmPickedHelperStructs::MassTrace MassTrace;
  typedef FeatureFinderAlgorithmPickedHelperStructs::MassTraces MassTraces;

  typedef std::multimap<double, PeptideIdentification*> RTMap;
  typedef std::map<Int, std::pair<RTMap, RTMap> > ChargeMap;
  typedef std::map<AASequence, ChargeMap> PeptideMap;
  typedef std::map<String, std::pair<RTMap, RTMap> > PeptideRefRTMap;

  PeptideMap peptide_map_;

  Size n_internal_peps_; 
  Size n_external_peps_; 

  Size batch_size_; 
  double rt_window_; 
  double mz_window_; 
  bool mz_window_ppm_; 

  double mapping_tolerance_; 

  double isotope_pmin_; 
  Size n_isotopes_; 

  double rt_quantile_;

  double peak_width_;
  double min_peak_width_;
  double signal_to_noise_;

  String elution_model_;

  // SVM related parameters
  double svm_min_prob_;
  StringList svm_predictor_names_;
  String svm_xval_out_;
  double svm_quality_cutoff;
  Size svm_n_parts_; 
  Size svm_n_samples_; 

  // output file (before filtering)
  String candidates_out_;

  Size debug_level_;

  void updateMembers_() override;

  struct RTRegion
  {
    double start, end;
    ChargeMap ids; 
  };

  struct FeatureFilterQuality
  {
    bool operator()(const Feature& feature)
    {
      return feature.getOverallQuality() == 0.0;
    }
  } feature_filter_quality_;

  struct FeatureFilterPeptides
  {
    bool operator()(const Feature& feature)
    {
      return feature.getPeptideIdentifications().empty();
    }
  } feature_filter_peptides_;

  struct PeptideCompare
  {
    bool operator()(const PeptideIdentification& p1,
                    const PeptideIdentification& p2)
    {
      const String& seq1 = p1.getHits()[0].getSequence().toString();
      const String& seq2 = p2.getHits()[0].getSequence().toString();
      if (seq1 == seq2)
      {
        Int charge1 = p1.getHits()[0].getCharge();
        Int charge2 = p2.getHits()[0].getCharge();
        if (charge1 == charge2)
        {
          return p1.getRT() < p2.getRT();
        }
        return charge1 < charge2;
      }
      return seq1 < seq2;
    }
  } peptide_compare_;

  struct FeatureCompare
  {
    bool operator()(const Feature& f1, const Feature& f2)
    {
      const String& ref1 = f1.getMetaValue("PeptideRef");
      const String& ref2 = f2.getMetaValue("PeptideRef");
      if (ref1 == ref2)
      {
        return f1.getRT() < f2.getRT();
      }
      return ref1 < ref2;
    }
  } feature_compare_;

  PeakMap ms_data_; 
  PeakMap chrom_data_; 
  TargetedExperiment library_; 

  bool quantify_decoys_;
  bool use_psm_cutoff_;
  double psm_score_cutoff_;
  std::vector<PeptideIdentification> unassignedIDs_;

  const double seed_rt_window_ = 60.0; 

  std::map<double, std::pair<Size, Size> > svm_probs_internal_;
  std::multiset<double> svm_probs_external_;
  Size n_internal_features_; 
  Size n_external_features_; 
  TransformationDescription trafo_external_; 
  std::map<String, double> isotope_probs_; 
  MRMFeatureFinderScoring feat_finder_; 

  ProgressLogger prog_log_;

  void generateTransitions_(const String& peptide_id, double mz, Int charge,
                            const IsotopeDistribution& iso_dist);

  void addPeptideRT_(TargetedExperiment::Peptide& peptide, double rt) const;

  void getRTRegions_(ChargeMap& peptide_data, std::vector<RTRegion>& rt_regions, bool clear_IDs = true) const;

  void annotateFeaturesFinalizeAssay_(
    FeatureMap& features,
    std::map<Size, std::vector<PeptideIdentification*> >& feat_ids,
    RTMap& rt_internal);

  void annotateFeatures_(FeatureMap& features, PeptideRefRTMap& ref_rt_map);

  void ensureConvexHulls_(Feature& feature) const;

  void postProcess_(FeatureMap& features, bool with_external_ids);

  void statistics_(const FeatureMap& features) const;

  void createAssayLibrary_(const PeptideMap::iterator& begin, const PeptideMap::iterator& end, PeptideRefRTMap& ref_rt_map, bool clear_IDs = true);

  void addPeptideToMap_(PeptideIdentification& peptide,
    PeptideMap& peptide_map,
    bool external = false);

  void checkNumObservations_(Size n_pos, Size n_neg, const String& note = "") const;

  void getUnbiasedSample_(const std::multimap<double, std::pair<Size, bool> >& valid_obs,
                          std::map<Size, double>& training_labels);

  void getRandomSample_(std::map<Size, double>& training_labels) const;

  void classifyFeatures_(FeatureMap& features);

  void filterFeaturesFinalizeAssay_(Feature& best_feature, double best_quality,
                                    const double quality_cutoff);

  void filterFeatures_(FeatureMap& features, bool classified);

  void calculateFDR_(FeatureMap& features);

  template <typename It>
  std::vector<std::pair<It,It>>
  chunk_(It range_from, It range_to, const std::ptrdiff_t batch_size)
  {
    /* Aliases, to make the rest of the code more readable. */
    using std::vector;
    using std::pair;
    using std::make_pair;
    using std::distance;
    using diff_t = std::ptrdiff_t;

    /* Total item number and batch_size size. */
    const diff_t total {distance(range_from, range_to)};
    const diff_t num {total / batch_size};

    vector<pair<It,It>> chunks(num);

    It batch_end {range_from};

    /* Use the 'generate' algorithm to create batches. */
    std::generate(begin(chunks), end(chunks), [&batch_end, batch_size]()
    {
      It batch_start {batch_end };

      std::advance(batch_end, batch_size);
      return make_pair(batch_start, batch_end);
    });

    /* The last batch_size's end must always be 'range_to'. */
    if (chunks.empty())
    {
      chunks.emplace_back(range_from, range_to);
    }
    else
    {
      chunks.back().second = range_to;
    }

    return chunks;
  }
};

} // namespace OpenMS

#endif