OpenMS
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
MapAlignmentAlgorithmIdentification.h
Go to the documentation of this file.
1 // Copyright (c) 2002-present, OpenMS Inc. -- EKU Tuebingen, ETH Zurich, and FU Berlin
2 // SPDX-License-Identifier: BSD-3-Clause
3 //
4 // --------------------------------------------------------------------------
5 // $Maintainer: Hendrik Weisser $
6 // $Authors: Eva Lange, Clemens Groepl, Hendrik Weisser $
7 // --------------------------------------------------------------------------
8 
9 #pragma once
10 
22 
23 #include <cmath> // for "abs"
24 #include <limits> // for "max"
25 #include <map>
26 
27 namespace OpenMS
28 {
29  /* Concept for FeatureMap or ConsensusMap*/
30  template <typename MapType>
31  concept IsFCMap = std::same_as<MapType, OpenMS::FeatureMap> || std::same_as<MapType, OpenMS::ConsensusMap>;
32 
33  class AnnotatedMSRun;
34 
54  public DefaultParamHandler,
55  public ProgressLogger
56  {
57 public:
60 
63 
64  // Set a reference for the alignment
65  template <typename DataType> void setReference(const DataType& data)
66  {
67  reference_.clear();
68  if (data.empty()) return; // empty input resets the reference
69  SeqToList rt_data;
70  // set these here because "checkParameters_" may not have been called yet:
71  use_feature_rt_ = param_.getValue("use_feature_rt").toBool();
72  score_cutoff_ = param_.getValue("score_cutoff").toBool();
73  score_type_ = (std::string)param_.getValue("score_type");
74  bool sorted = getRetentionTimes_(data, rt_data);
75  computeMedians_(rt_data, reference_, sorted);
76 
77  if (reference_.empty())
78  {
79  throw Exception::MissingInformation(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "Could not extract retention time information from the reference file");
80  }
81  }
82 
92  template <typename DataType>
93  void align(const std::vector<DataType>& data,
94  std::vector<TransformationDescription>& transformations,
95  Int reference_index = -1)
96  {
97  checkParameters_(data.size());
98  startProgress(0, 3, "aligning maps");
99 
100  reference_index_ = reference_index;
101  // is reference one of the input files?
102  bool use_internal_reference = (reference_index >= 0);
103  if (use_internal_reference)
104  {
105  if (reference_index >= Int(data.size()))
106  {
107  throw Exception::IndexOverflow(__FILE__, __LINE__,
108  OPENMS_PRETTY_FUNCTION,
109  reference_index, data.size());
110  }
111  setReference(data[reference_index]);
112  }
113 
114  // one set of RT data for each input map, except reference (if any):
115  std::vector<SeqToList> rt_data(data.size() - use_internal_reference);
116  bool all_sorted = true;
117  for (Size i = 0, j = 0; i < data.size(); ++i)
118  {
119  if ((reference_index >= 0) && (i == Size(reference_index)))
120  {
121  continue; // skip reference map, if any
122  }
123  all_sorted &= getRetentionTimes_(data[i], rt_data[j++]);
124  }
125  setProgress(1);
126 
127  computeTransformations_(rt_data, transformations, all_sorted);
128  setProgress(2);
129 
130  setProgress(3);
131  endProgress();
132  }
133 
134 protected:
135 
137  typedef std::map<String, DoubleList> SeqToList;
138 
140  typedef std::map<String, double> SeqToValue;
141 
144 
147 
150 
152  bool use_feature_rt_{};
153 
155  bool use_adducts_{};
156 
158  double min_score_;
159 
161  bool score_cutoff_{};
162 
165 
167  bool (*better_) (double, double) = [](double, double) {return true;};
168 
178  void computeMedians_(SeqToList& rt_data, SeqToValue& medians,
179  bool sorted = false);
180 
190  SeqToList& rt_data);
191 
200  // "id_data" can't be "const" here or template resolution will fail
201  bool getRetentionTimes_(const IdentificationData& id_data, SeqToList& rt_data);
202 
218  bool getRetentionTimes_(const IsFCMap auto& features, SeqToList& rt_data)
219  {
220  if (!score_cutoff_)
221  {
222  better_ = [](double, double)
223  {return true;};
224  }
225  else if (features[0].getPeptideIdentifications()[0].isHigherScoreBetter())
226  {
227  better_ = [](double a, double b)
228  { return a >= b; };
229  }
230  else
231  {
232  better_ = [](double a, double b)
233  { return a <= b; };
234  }
235 
236  for (auto feat_it = features.cbegin(); feat_it != features.cend(); ++feat_it)
237  {
238  if (use_feature_rt_)
239  {
240  // find the peptide ID closest in RT to the feature centroid:
241  String sequence;
242  double rt_distance = std::numeric_limits<double>::max();
243  bool any_hit = false;
245  feat_it->getPeptideIdentifications().begin(); pep_it !=
246  feat_it->getPeptideIdentifications().end(); ++pep_it)
247  {
248  if (!pep_it->getHits().empty())
249  {
250  any_hit = true;
251  double current_distance = fabs(pep_it->getRT() -
252  feat_it->getRT());
253  if (current_distance < rt_distance)
254  {
255  const PeptideHit* best_hit = getBestScoringHit(pep_it->getHits(), pep_it->isHigherScoreBetter());
256  if (best_hit && better_(best_hit->getScore(), min_score_))
257  {
258  sequence = best_hit->getSequence().toString();
259  rt_distance = current_distance;
260  }
261  }
262  }
263  }
264 
265  if (any_hit) rt_data[sequence].push_back(feat_it->getRT());
266  }
267  else
268  {
269  getRetentionTimes_(feat_it->getPeptideIdentifications(), rt_data);
270  }
271  }
272 
273  if (!use_feature_rt_ &&
274  param_.getValue("use_unassigned_peptides").toBool())
275  {
276  getRetentionTimes_(features.getUnassignedPeptideIdentifications(),
277  rt_data);
278  }
279 
280  // remove duplicates (can occur if a peptide ID was assigned to several
281  // features due to overlap or annotation tolerance):
282  for (SeqToList::iterator rt_it = rt_data.begin(); rt_it != rt_data.end();
283  ++rt_it)
284  {
285  DoubleList& rt_values = rt_it->second;
286  sort(rt_values.begin(), rt_values.end());
287  DoubleList::iterator it = unique(rt_values.begin(), rt_values.end());
288  rt_values.resize(it - rt_values.begin());
289  }
290  return true; // RTs were already sorted for duplicate detection
291  }
292 
300  void computeTransformations_(std::vector<SeqToList>& rt_data,
301  std::vector<TransformationDescription>&
302  transforms, bool sorted = false);
303 
311  void checkParameters_(const Size runs);
312 
319 
326 
335  const PeptideHit* getBestScoringHit(const std::vector<PeptideHit>& hits, const bool is_higher_score_better);
336 
337 private:
338 
341 
344 
345  };
346 
347 } // namespace OpenMS
String toString() const
returns the peptide as string with modifications embedded in brackets
Class for storing MS run data with peptide and protein identifications.
Definition: AnnotatedMSRun.h:36
A base class for all classes handling default parameters.
Definition: DefaultParamHandler.h:66
Int overflow exception.
Definition: Exception.h:211
Not all required information provided.
Definition: Exception.h:155
typename VecMember::const_iterator const_iterator
Definition: ExposedVector.h:69
Definition: IdentificationData.h:87
A map alignment algorithm based on peptide identifications from MS2 spectra.
Definition: MapAlignmentAlgorithmIdentification.h:56
void computeTransformations_(std::vector< SeqToList > &rt_data, std::vector< TransformationDescription > &transforms, bool sorted=false)
Compute retention time transformations from RT data grouped by peptide sequence.
bool getRetentionTimes_(const PeptideIdentificationList &peptides, SeqToList &rt_data)
Collect retention time data from peptide IDs.
void setReference(const DataType &data)
Definition: MapAlignmentAlgorithmIdentification.h:65
bool getRetentionTimes_(const IdentificationData &id_data, SeqToList &rt_data)
Collect retention time data from spectrum matches.
~MapAlignmentAlgorithmIdentification() override
Destructor.
std::map< String, double > SeqToValue
Type to store one representative retention time per peptide sequence.
Definition: MapAlignmentAlgorithmIdentification.h:140
void checkParameters_(const Size runs)
Check that parameter values are valid.
void getReference_()
Get reference retention times.
bool getRetentionTimes_(const IsFCMap auto &features, SeqToList &rt_data)
Collect retention time data from peptide IDs contained in feature maps or consensus maps.
Definition: MapAlignmentAlgorithmIdentification.h:218
String score_type_
Score type to use for filtering.
Definition: MapAlignmentAlgorithmIdentification.h:164
Int reference_index_
Index of input file to use as reference (if any)
Definition: MapAlignmentAlgorithmIdentification.h:143
const PeptideHit * getBestScoringHit(const std::vector< PeptideHit > &hits, const bool is_higher_score_better)
Get the best-scoring PeptideHit from a list of hits.
SeqToValue reference_
Reference retention times (per peptide sequence)
Definition: MapAlignmentAlgorithmIdentification.h:146
double min_score_
Minimum score to reach for a peptide to be considered.
Definition: MapAlignmentAlgorithmIdentification.h:158
Size min_run_occur_
Minimum number of runs a peptide must occur in.
Definition: MapAlignmentAlgorithmIdentification.h:149
std::map< String, DoubleList > SeqToList
Type to store retention times given for individual peptide sequences.
Definition: MapAlignmentAlgorithmIdentification.h:137
MapAlignmentAlgorithmIdentification & operator=(const MapAlignmentAlgorithmIdentification &)
Assignment operator intentionally not implemented -> private.
MapAlignmentAlgorithmIdentification()
Default constructor.
IdentificationData::ScoreTypeRef handleIdDataScoreType_(const IdentificationData &id_data)
Helper function to find/define the score type for processing IdentificationData.
void align(const std::vector< DataType > &data, std::vector< TransformationDescription > &transformations, Int reference_index=-1)
Align feature maps, consensus maps, or peptide identifications.
Definition: MapAlignmentAlgorithmIdentification.h:93
void computeMedians_(SeqToList &rt_data, SeqToValue &medians, bool sorted=false)
Compute the median retention time for each peptide sequence.
MapAlignmentAlgorithmIdentification(const MapAlignmentAlgorithmIdentification &)
Copy constructor intentionally not implemented -> private.
Represents a single spectrum match (candidate) for a specific tandem mass spectrum (MS/MS).
Definition: PeptideHit.h:50
double getScore() const
returns the PSM score
const AASequence & getSequence() const
returns the peptide sequence
Container for peptide identifications from multiple spectra.
Definition: PeptideIdentificationList.h:66
Base class for all classes that want to report their progress.
Definition: ProgressLogger.h:27
A more convenient string class.
Definition: String.h:34
int Int
Signed integer type.
Definition: Types.h:72
size_t Size
Size type e.g. used as variable which can hold result of size()
Definition: Types.h:97
std::vector< double > DoubleList
Vector of double precision real types.
Definition: ListUtils.h:36
Main OpenMS namespace.
Definition: openswathalgo/include/OpenMS/OPENSWATHALGO/DATAACCESS/ISpectrumAccess.h:19
concept IsFCMap
Definition: MapAlignmentAlgorithmIdentification.h:31
Wrapper that adds operator< to iterators, so they can be used as (part of) keys in maps/sets or multi...
Definition: MetaData.h:20