35 #ifndef OPENMS_ANALYSIS_TARGETED_OFFLINEPRECURSORIONSELECTION_H
36 #define OPENMS_ANALYSIS_TARGETED_OFFLINEPRECURSORIONSELECTION_H
48 class PeptideIdentification;
49 class ProteinIdentification;
80 template <
typename InputPeakType>
81 void makePrecursorSelectionForKnownLCMSMap(
const FeatureMap& features,
84 std::set<Int>& charges_set,
94 template <
typename InputPeakType>
97 std::vector<std::vector<std::pair<Size, Size> > >& indices);
104 std::cout <<
" LPSolver set to " << solver_ << std::endl;
116 template <
typename InputPeakType>
117 void calculateXICs_(
const FeatureMap& features,
118 const std::vector<std::vector<std::pair<Size, Size> > >& mass_ranges,
120 const std::set<Int>& charges_set,
121 std::vector<std::vector<std::pair<Size, double> > >& xics);
126 template <
typename InputPeakType>
127 void checkMassRanges_(std::vector<std::vector<std::pair<Size, Size> > >& mass_ranges,
130 template <
typename T>
131 void updateExclusionList_(std::vector<std::pair<T, Size> >& exclusion_list);
138 template <
typename InputPeakType>
141 bool enclose_hit =
false;
143 for (
Size i = 0; i < hulls.size(); ++i)
145 if (hulls[i].getBoundingBox().encloses(rt, mz))
154 template <
typename InputPeakType>
157 std::vector<std::vector<std::pair<Size, Size> > >& indices)
159 if (experiment.
empty())
161 for (
Size f = 0; f < features.size(); ++f)
163 std::vector<std::pair<Size, Size> > vec;
165 for (
Size rt = 0; rt < experiment.
size(); ++rt)
168 if (!enclosesBoundingBox<InputPeakType>(features[f], experiment[rt].getRT(), features[f].getMZ()))
171 std::pair<Size, Size> start;
172 std::pair<Size, Size> end;
173 bool start_found =
false;
174 bool end_found =
false;
177 if (mz_iter == experiment[rt].end())
180 while (enclosesBoundingBox<InputPeakType>(features[f], experiment[rt].getRT(), mz_iter->getMZ()))
184 start.second = distance(experiment[rt].begin(), mz_iter);
185 if (mz_iter == experiment[rt].begin())
190 end.second = start.second;
193 while (mz_end != experiment[rt].end() && enclosesBoundingBox<InputPeakType>(features[f], experiment[rt].getRT(), mz_end->getMZ()))
197 end.second = distance(experiment[rt].begin(), mz_end);
200 if (start_found && end_found)
202 vec.push_back(start);
206 else if (start_found || end_found)
208 std::cout <<
"start " << start_found <<
" end " << end_found << std::endl;
209 std::cout <<
"feature: " << f <<
" rt: " << rt << std::endl;
216 std::cout << vec.size() <<
" / 2 scans" << std::endl;
217 for (
Size i = 0; i < vec.size(); i += 2)
219 std::cout <<
"Feature " << f <<
" RT : " << vec[i].first
220 <<
" MZ : " << experiment[vec[i].first][vec[i].second].getMZ() <<
" "
221 << experiment[vec[i + 1].first][vec[i + 1].second].getMZ() << std::endl;
228 std::cout <<
"According to the convex hulls no mass traces found for this feature->estimate!"
229 << features[f].getRT() <<
" " << features[f].getMZ() <<
" " << features[f].getCharge() << std::endl;
233 if (spec_iter == experiment.
end())
236 double dist1 = fabs(spec_iter->getRT() - features[f].getRT());
237 double dist2 = std::numeric_limits<double>::max();
238 double dist3 = std::numeric_limits<double>::max();
239 if ((spec_iter + 1) != experiment.
end())
241 dist2 = fabs((spec_iter + 1)->getRT() - features[f].getRT());
243 if (spec_iter != experiment.
begin())
245 dist3 = fabs((spec_iter - 1)->getRT() - features[f].getRT());
247 if (dist3 <= dist1 && dist3 <= dist2)
251 else if (dist2 <= dist3 && dist2 <= dist1)
255 std::pair<Size, Size> start;
256 std::pair<Size, Size> end;
257 start.first = distance(experiment.
begin(), spec_iter);
258 end.first = start.first;
261 if (spec_iter->
begin() == spec_iter->
end())
263 indices.push_back(vec);
266 if (mz_iter == spec_iter->
end() || (mz_iter->getMZ() > features[f].getMZ() && mz_iter != spec_iter->
begin()))
268 while (mz_iter != spec_iter->
begin())
270 if (fabs((mz_iter - 1)->getMZ() - features[f].getMZ()) < 0.5)
275 start.second = distance(spec_iter->
begin(), mz_iter);
278 std::cout << features[f].getMZ() <<
" Start: " << experiment[start.first].getRT() <<
" " << experiment[start.first][start.second].getMZ();
280 Int charge = features[f].getCharge();
283 while (mz_end + 1 != spec_iter->
end())
285 if (fabs((mz_end + 1)->getMZ() - features[f].getMZ()) < 3.0 / (
double)charge)
290 end.second = distance(spec_iter->
begin(), mz_end);
292 std::cout <<
"\tEnd: " << experiment[end.first].getRT() <<
" " << experiment[end.first][end.second].getMZ() << std::endl;
294 vec.push_back(start);
298 indices.push_back(vec);
305 template <
typename InputPeakType>
307 const std::vector<std::vector<std::pair<Size, Size> > >& mass_ranges,
309 const std::set<Int>& charges_set,
310 std::vector<std::vector<std::pair<Size, double> > >& xics)
313 xics.resize(experiment.
size());
315 for (
Size f = 0; f < mass_ranges.size(); ++f)
318 if (charges_set.count(features[f].getCharge()) < 1)
323 for (
Size s = 0; s < mass_ranges[f].size(); s += 2)
327 for (
Size j = mass_ranges[f][s].second; j <= mass_ranges[f][s + 1].second; ++j)
329 weight += experiment[mass_ranges[f][s].first][j].getIntensity();
332 xics[mass_ranges[f][s].first].push_back(std::make_pair(f, weight));
336 for (
Size s = 0; s < xics.size(); ++s)
342 template <
typename InputPeakType>
346 std::set<Int>& charges_set,
354 std::vector<std::vector<std::pair<Size, Size> > > indices;
357 if (experiment.
size() > 1)
359 rt_dist = experiment[1].getRT() - experiment[0].getRT();
368 std::vector<IndexTriple> variable_indices;
369 std::vector<int> solution_indices;
371 indices, charges_set,
377 std::cout <<
"best_solution " << std::endl;
381 for (
Size i = 0; i < solution_indices.size(); ++i)
383 Size feature_index = variable_indices[solution_indices[i]].feature;
384 Size feature_scan_idx = variable_indices[solution_indices[i]].scan;
388 std::vector<Precursor> pcs;
390 p.
setMZ(features[feature_index].getMZ());
391 p.
setCharge(features[feature_index].getCharge());
394 ms2_spec.
setRT(scan->getRT() + rt_dist / 2.0);
399 ms2_spec.
setMetaValue(
"parent_feature_ids", ListUtils::create<Int>(
String(feature_index)));
401 std::cout <<
" MS2 spectra generated at: " << scan->getRT() <<
" x " << p.
getMZ() <<
"\n";
405 std::cout << solution_indices.size() <<
" out of " << features.size()
406 <<
" precursors are in best solution.\n";
412 std::cout <<
"scan based precursor selection" << std::endl;
417 std::vector<DoubleList> feature_elution_bounds;
418 std::vector<DoubleList> elution_profile_intensities;
419 std::vector<DoubleList> isotope_intensities;
421 bool meta_values_present =
false;
423 if (!features.empty() &&
424 features[0].metaValueExists(
"elution_profile_bounds") &&
425 features[0].metaValueExists(
"elution_profile_intensities") &&
426 features[0].metaValueExists(
"isotope_intensities"))
428 for (
Size feat = 0; feat < features.size(); ++feat)
430 feature_elution_bounds.push_back(features[feat].getMetaValue(
"elution_profile_bounds"));
431 elution_profile_intensities.push_back(features[feat].getMetaValue(
"elution_profile_intensities"));
432 isotope_intensities.push_back(features[feat].getMetaValue(
"isotope_intensities"));
434 meta_values_present =
true;
438 std::vector<std::vector<Size> > scan_features(experiment.
size());
440 for (
Size feat = 0; feat < features.size(); ++feat)
442 if (charges_set.count(features[feat].getCharge()))
444 Size lower_rt = features[feat].getConvexHull().getBoundingBox().minX();
445 Size upper_rt = features[feat].getConvexHull().getBoundingBox().maxX();
447 for (it = experiment.
RTBegin(lower_rt); it != experiment.
RTEnd(upper_rt); ++it)
449 scan_features[it - experiment.
begin()].push_back(feat);
454 bool dynamic_exclusion =
param_.
getValue(
"Exclusion:use_dynamic_exclusion") ==
"true" ?
true :
false;
456 ExclusionListType exclusion_list;
458 if (!dynamic_exclusion)
465 std::map<Size, typename OpenMS::DBoundingBox<2> > bounding_boxes_f;
467 for (Size feature_num = 0; feature_num < features.size(); ++feature_num)
469 if (charges_set.count(features[feature_num].getCharge()))
471 bounding_boxes_f.insert(std::make_pair(feature_num, features[feature_num].getConvexHull().getBoundingBox()));
472 const std::vector<ConvexHull2D> mass_traces = features[feature_num].getConvexHulls();
473 for (Size mass_trace_num = 0; mass_trace_num < mass_traces.size(); ++mass_trace_num)
478 bounding_boxes.insert(std::make_pair(std::make_pair(feature_num, mass_trace_num), tmp_bbox));
485 for (Size i = 0; i < experiment.
size(); ++i)
488 std::cout <<
"scan " << experiment[i].getRT() <<
":";
494 Size selected_peaks = 0, j = 0;
496 while (selected_peaks < max_spec && j < scan.size())
498 double peak_mz = scan[j].getMZ();
499 double peak_rt = scan.
getRT();
501 ExclusionListType::iterator it_low = exclusion_list.lower_bound(std::make_pair(peak_mz, peak_mz));
502 if (it_low != exclusion_list.end() && it_low->first.first <= peak_mz)
511 std::vector<Precursor> pcs;
512 std::set<std::pair<Size, Size> > selected_mt;
515 double local_mz = peak_mz;
517 for (Size scan_feat_id = 0; scan_feat_id < scan_features[i].size(); ++scan_feat_id)
519 Size feature_num = scan_features[i][scan_feat_id];
520 if (bounding_boxes_f[feature_num].encloses(peak_rt, local_mz))
523 double feature_intensity = 0;
524 for (Size mass_trace_num = 0; mass_trace_num < features[feature_num].getConvexHulls().size(); ++mass_trace_num)
526 if (bounding_boxes[std::make_pair(feature_num, mass_trace_num)].encloses(
DPosition<2>(peak_rt, local_mz)))
528 double elu_factor = 1.0, iso_factor = 1.0;
530 if (meta_values_present)
532 DoubleList xxx = elution_profile_intensities[feature_num];
533 DoubleList yyy = feature_elution_bounds[feature_num];
537 OPENMS_PRECONDITION(i - yyy[0] < xxx.size(),
"Tried to access invalid index for elution factor");
538 elu_factor = xxx[i - yyy[0]];
539 iso_factor = isotope_intensities[feature_num][mass_trace_num];
541 feature_intensity += features[feature_num].getIntensity() * iso_factor * elu_factor;
546 p.
setMZ(features[feature_num].getMZ());
547 p.
setCharge(features[feature_num].getCharge());
549 parent_feature_ids.push_back((
Int)feature_num);
558 ms2_spec.
setRT(experiment[i].getRT() + rt_dist / 2.0);
559 ms2_spec.
setMetaValue(
"parent_feature_ids", parent_feature_ids);
564 exclusion_list.insert(std::make_pair(std::make_pair(peak_mz - excl_window, peak_mz + excl_window), exclusion_specs + 1));
572 template <
typename InputPeakType>
576 std::vector<std::vector<std::pair<Size, Size> > > checked_mass_ranges;
578 checked_mass_ranges.reserve(mass_ranges.size());
579 for (
Size f = 0; f < mass_ranges.size(); ++f)
581 std::vector<std::pair<Size, Size> > checked_mass_ranges_f;
582 for (
Size s_idx = 0; s_idx < mass_ranges[f].size(); s_idx += 2)
584 Size s = mass_ranges[f][s_idx].first;
585 bool overlapping_features =
false;
589 const InputPeakType& peak_left_border = experiment[s][mass_ranges[f][s_idx].second];
590 const InputPeakType& peak_right_border = experiment[s][mass_ranges[f][s_idx + 1].second];
591 for (
Size fmr = 0; fmr < mass_ranges.size(); ++fmr)
595 for (
Size mr = 0; mr < mass_ranges[fmr].size(); mr += 2)
597 if (mass_ranges[fmr][mr].first == s)
599 const InputPeakType& tmp_peak_left = experiment[s][mass_ranges[fmr][mr].second];
600 const InputPeakType& tmp_peak_right = experiment[s][mass_ranges[fmr][mr + 1].second];
602 std::cout << tmp_peak_left.getMZ() <<
" < "
603 << peak_left_border.getMZ() - min_peak_distance <<
" && "
604 << tmp_peak_right.getMZ() <<
" < "
605 << peak_left_border.getMZ() - min_peak_distance <<
" ? "
606 << (tmp_peak_left.getMZ() < peak_left_border.getMZ() - min_peak_distance &&
607 tmp_peak_right.getMZ() < peak_left_border.getMZ() - min_peak_distance)
609 << tmp_peak_left.getMZ() <<
" > "
610 << peak_right_border.getMZ() + min_peak_distance <<
" && "
611 << tmp_peak_right.getMZ() <<
" > "
612 << peak_right_border.getMZ() + min_peak_distance <<
" ? "
613 << (tmp_peak_left.getMZ() > peak_right_border.getMZ() + min_peak_distance &&
614 tmp_peak_right.getMZ() > peak_right_border.getMZ() + min_peak_distance)
619 if (!((tmp_peak_left.getMZ() < peak_left_border.getMZ() - min_peak_distance &&
620 tmp_peak_right.getMZ() < peak_left_border.getMZ() - min_peak_distance) ||
621 (tmp_peak_left.getMZ() > peak_right_border.getMZ() + min_peak_distance &&
622 tmp_peak_right.getMZ() > peak_right_border.getMZ() + min_peak_distance)))
625 std::cout <<
"found overlapping peak" << std::endl;
627 overlapping_features =
true;
633 if (!overlapping_features)
636 std::cout <<
"feature in spec ok" << mass_ranges[f][s_idx].second <<
" in spec "
637 << mass_ranges[f][s_idx].first << std::endl;
639 checked_mass_ranges_f.insert(checked_mass_ranges_f.end(),
640 mass_ranges[f].begin() + s_idx,
641 mass_ranges[f].begin() + s_idx + 2);
644 checked_mass_ranges.push_back(checked_mass_ranges_f);
646 mass_ranges.swap(checked_mass_ranges);
649 template <
typename T>
652 for (
Size i = 0; i < exclusion_list.size(); ++i)
654 if (exclusion_list[i].second > 0)
655 --exclusion_list[i].second;
658 typename std::vector<std::pair<T, Size> >::iterator iter = exclusion_list.begin();
659 while (iter != exclusion_list.end() && iter->second != 0)
661 exclusion_list.erase(iter, exclusion_list.end());
668 it = exclusion_list.begin();
670 while (it != exclusion_list.end())
672 if ((it->second--) == 1)
674 exclusion_list.erase(it++);
685 #endif // OPENMS_ANALYSIS_ID_OFFLINEPRECURSORIONSELECTION_H
CoordinateType maxY() const
Accessor for max_ coordinate maximum.
Definition: DIntervalBase.h:259
A more convenient string class.
Definition: String.h:57
Precursor meta information.
Definition: Precursor.h:56
void setMinY(CoordinateType const c)
Mutator for max_ coordinate of the smaller point.
Definition: DIntervalBase.h:272
Size size() const
Definition: MSExperiment.h:117
std::vector< double > DoubleList
Vector of double precision real types.
Definition: ListUtils.h:66
Param param_
Container for current parameters.
Definition: DefaultParamHandler.h:150
SOLVER
Definition: LPWrapper.h:129
A container for features.
Definition: FeatureMap.h:93
#define OPENMS_PRECONDITION(condition, message)
Precondition macro.
Definition: openms/include/OpenMS/CONCEPT/Macros.h:107
void setMaxY(CoordinateType const c)
Mutator for max_ coordinate of the larger point.
Definition: DIntervalBase.h:286
LPWrapper::SOLVER getLPSolver()
Definition: OfflinePrecursorIonSelection.h:107
CoordinateType getMZ() const
Non-mutable access to m/z.
Definition: Peak1D.h:114
ContainerType::const_iterator ConstIterator
Non-mutable iterator.
Definition: MSSpectrum.h:123
Iterator begin()
Definition: MSExperiment.h:147
bool enclosesBoundingBox(const Feature &f, typename MSExperiment< InputPeakType >::CoordinateType rt, typename MSExperiment< InputPeakType >::CoordinateType mz)
Definition: OfflinePrecursorIonSelection.h:139
std::vector< Int > IntList
Vector of signed integers.
Definition: ListUtils.h:59
Main OpenMS namespace.
Definition: FeatureDeconvolution.h:47
void setMZ(CoordinateType mz)
Mutable access to m/z.
Definition: Peak1D.h:120
void setIntensity(IntensityType intensity)
Mutable access to the data point intensity (height)
Definition: Peak1D.h:111
ConstIterator RTEnd(CoordinateType rt) const
Fast search for spectrum range end (returns the past-the-end iterator)
Definition: MSExperiment.h:388
Class for comparison of std::pair using second ONLY e.g. for use with std::sort.
Definition: ComparatorUtils.h:368
void makePrecursorSelectionForKnownLCMSMap(const FeatureMap &features, const MSExperiment< InputPeakType > &experiment, MSExperiment< InputPeakType > &ms2, std::set< Int > &charges_set, bool feature_based)
Makes the precursor selection for a given feature map, either feature or scan based.
Definition: OfflinePrecursorIonSelection.h:343
Iterator end()
Definition: MSExperiment.h:157
ConstIterator RTBegin(CoordinateType rt) const
Fast search for spectrum range begin.
Definition: MSExperiment.h:374
CoordinateType minY() const
Accessor for max_ coordinate minimum.
Definition: DIntervalBase.h:247
void getMassRanges(const FeatureMap &features, const MSExperiment< InputPeakType > &experiment, std::vector< std::vector< std::pair< Size, Size > > > &indices)
Calculates the mass ranges for each feature and stores them as indices of the raw data...
Definition: OfflinePrecursorIonSelection.h:155
const DataValue & getValue(const String &key) const
Returns a value of a parameter.
void calculateXICs_(const FeatureMap &features, const std::vector< std::vector< std::pair< Size, Size > > > &mass_ranges, const MSExperiment< InputPeakType > &experiment, const std::set< Int > &charges_set, std::vector< std::vector< std::pair< Size, double > > > &xics)
Calculate the sum of intensities of relevant features for each scan separately.
Definition: OfflinePrecursorIonSelection.h:306
void updateExclusionList_(std::vector< std::pair< T, Size > > &exclusion_list)
Definition: OfflinePrecursorIonSelection.h:650
const std::vector< ConvexHull2D > & getConvexHulls() const
Non-mutable access to the convex hulls.
void sortByIntensity(bool reverse=false)
Lexicographically sorts the peaks by their intensity.
Definition: MSSpectrum.h:342
LPWrapper::SOLVER solver_
Definition: OfflinePrecursorIonSelection.h:135
double getRT() const
Definition: MSSpectrum.h:243
void setMSLevel(UInt ms_level)
Sets the MS level.
Definition: MSSpectrum.h:265
An LC-MS feature.
Definition: Feature.h:70
void setLPSolver(LPWrapper::SOLVER solver)
Definition: OfflinePrecursorIonSelection.h:101
void addSpectrum(const MSSpectrum< PeakT > &spectrum)
adds a spectra to the list
Definition: MSExperiment.h:758
void setRT(double rt)
Sets the absolute retention time (is seconds)
Definition: MSSpectrum.h:249
Invalid UInt exception.
Definition: Exception.h:304
void checkMassRanges_(std::vector< std::vector< std::pair< Size, Size > > > &mass_ranges, const MSExperiment< InputPeakType > &experiment)
Eliminates overlapping peaks.
Definition: OfflinePrecursorIonSelection.h:573
In-Memory representation of a mass spectrometry experiment.
Definition: MSExperiment.h:69
std::vector< SpectrumType >::const_iterator ConstIterator
Non-mutable iterator.
Definition: MSExperiment.h:103
Implements different algorithms for precursor ion selection.
Definition: OfflinePrecursorIonSelection.h:62
void setPrecursors(const std::vector< Precursor > &precursors)
sets the precursors
size_t Size
Size type e.g. used as variable which can hold result of size()
Definition: Types.h:121
Class for comparison of std::pair using second ONLY e.g. for use with std::sort.
Definition: ComparatorUtils.h:340
bool empty() const
Definition: MSExperiment.h:127
void setCharge(Int charge)
Mutable access to the charge.
A base class for all classes handling default parameters.
Definition: DefaultParamHandler.h:92
PSLPFormulation::IndexTriple IndexTriple
Definition: OfflinePrecursorIonSelection.h:66
int Int
Signed integer type.
Definition: Types.h:96