35 #ifndef OPENMS_ANALYSIS_TARGETED_OFFLINEPRECURSORIONSELECTION_H
36 #define OPENMS_ANALYSIS_TARGETED_OFFLINEPRECURSORIONSELECTION_H
47 class PeptideIdentification;
48 class ProteinIdentification;
79 template <
typename InputPeakType>
80 void makePrecursorSelectionForKnownLCMSMap(
const FeatureMap<>& features,
83 std::set<Int>& charges_set,
93 template <
typename InputPeakType>
96 std::vector<std::vector<std::pair<Size, Size> > >& indices);
103 std::cout <<
" LPSolver set to " << solver_ << std::endl;
115 template <
typename InputPeakType>
117 const std::vector<std::vector<std::pair<Size, Size> > >& mass_ranges,
119 const std::set<Int>& charges_set,
120 std::vector<std::vector<std::pair<Size, DoubleReal> > >& xics);
125 template <
typename InputPeakType>
126 void checkMassRanges_(std::vector<std::vector<std::pair<Size, Size> > >& mass_ranges,
129 template <
typename T>
130 void updateExclusionList_(std::vector<std::pair<T, Size> >& exclusion_list);
132 void updateExclusionList_(std::map<std::pair<DoubleReal, DoubleReal>,
Size,
PairComparatorSecondElement<std::pair<DoubleReal, DoubleReal> > >& exclusion_list);
137 template <
typename InputPeakType>
140 bool enclose_hit =
false;
142 for (
Size i = 0; i < hulls.size(); ++i)
144 if (hulls[i].getBoundingBox().encloses(rt, mz))
153 template <
typename InputPeakType>
156 std::vector<std::vector<std::pair<Size, Size> > >& indices)
158 if (experiment.
empty())
160 for (
Size f = 0; f < features.size(); ++f)
162 std::vector<std::pair<Size, Size> > vec;
164 for (
Size rt = 0; rt < experiment.
size(); ++rt)
167 if (!enclosesBoundingBox<InputPeakType>(features[f], experiment[rt].getRT(), features[f].getMZ()))
170 std::pair<Size, Size> start;
171 std::pair<Size, Size> end;
172 bool start_found =
false;
173 bool end_found =
false;
176 if (mz_iter == experiment[rt].end())
179 while (enclosesBoundingBox<InputPeakType>(features[f], experiment[rt].getRT(), mz_iter->getMZ()))
183 start.second = distance(experiment[rt].begin(), mz_iter);
184 if (mz_iter == experiment[rt].begin())
189 end.second = start.second;
192 while (mz_end != experiment[rt].end() && enclosesBoundingBox<InputPeakType>(features[f], experiment[rt].getRT(), mz_end->getMZ()))
196 end.second = distance(experiment[rt].begin(), mz_end);
199 if (start_found && end_found)
201 vec.push_back(start);
205 else if (start_found || end_found)
207 std::cout <<
"start " << start_found <<
" end " << end_found << std::endl;
208 std::cout <<
"feature: " << f <<
" rt: " << rt << std::endl;
215 std::cout << vec.size() <<
" / 2 scans" << std::endl;
216 for (
Size i = 0; i < vec.size(); i += 2)
218 std::cout <<
"Feature " << f <<
" RT : " << vec[i].first
219 <<
" MZ : " << experiment[vec[i].first][vec[i].second].getMZ() <<
" "
220 << experiment[vec[i + 1].first][vec[i + 1].second].getMZ() << std::endl;
227 std::cout <<
"According to the convex hulls no mass traces found for this feature->estimate!"
228 << features[f].getRT() <<
" " << features[f].getMZ() <<
" " << features[f].getCharge() << std::endl;
232 if (spec_iter == experiment.
end())
235 DoubleReal dist1 = fabs(spec_iter->getRT() - features[f].getRT());
236 DoubleReal dist2 = std::numeric_limits<DoubleReal>::max();
237 DoubleReal dist3 = std::numeric_limits<DoubleReal>::max();
238 if ((spec_iter + 1) != experiment.
end())
240 dist2 = fabs((spec_iter + 1)->getRT() - features[f].getRT());
242 if (spec_iter != experiment.
begin())
244 dist3 = fabs((spec_iter - 1)->getRT() - features[f].getRT());
246 if (dist3 <= dist1 && dist3 <= dist2)
250 else if (dist2 <= dist3 && dist2 <= dist1)
254 std::pair<Size, Size> start;
255 std::pair<Size, Size> end;
256 start.first = distance(experiment.
begin(), spec_iter);
257 end.first = start.first;
260 if (spec_iter->
begin() == spec_iter->
end())
262 indices.push_back(vec);
265 if (mz_iter == spec_iter->
end() || (mz_iter->getMZ() > features[f].getMZ() && mz_iter != spec_iter->
begin()))
267 while (mz_iter != spec_iter->
begin())
269 if (fabs((mz_iter - 1)->getMZ() - features[f].getMZ()) < 0.5)
274 start.second = distance(spec_iter->
begin(), mz_iter);
277 std::cout << features[f].getMZ() <<
" Start: " << experiment[start.first].getRT() <<
" " << experiment[start.first][start.second].getMZ();
279 Int charge = features[f].getCharge();
282 while (mz_end + 1 != spec_iter->
end())
284 if (fabs((mz_end + 1)->getMZ() - features[f].getMZ()) < 3.0 / (
DoubleReal)charge)
289 end.second = distance(spec_iter->
begin(), mz_end);
291 std::cout <<
"\tEnd: " << experiment[end.first].getRT() <<
" " << experiment[end.first][end.second].getMZ() << std::endl;
293 vec.push_back(start);
297 indices.push_back(vec);
304 template <
typename InputPeakType>
306 const std::vector<std::vector<std::pair<Size, Size> > >& mass_ranges,
308 const std::set<Int>& charges_set,
309 std::vector<std::vector<std::pair<Size, DoubleReal> > >& xics)
312 xics.resize(experiment.
size());
314 for (
Size f = 0; f < mass_ranges.size(); ++f)
317 if (charges_set.count(features[f].getCharge()) < 1)
322 for (
Size s = 0; s < mass_ranges[f].size(); s += 2)
326 for (
Size j = mass_ranges[f][s].second; j <= mass_ranges[f][s + 1].second; ++j)
328 weight += experiment[mass_ranges[f][s].first][j].getIntensity();
331 xics[mass_ranges[f][s].first].push_back(std::make_pair(f, weight));
335 for (
Size s = 0; s < xics.size(); ++s)
341 template <
typename InputPeakType>
345 std::set<Int>& charges_set,
353 std::vector<std::vector<std::pair<Size, Size> > > indices;
356 if (experiment.
size() > 1)
358 rt_dist = experiment[1].getRT() - experiment[0].getRT();
367 std::vector<IndexTriple> variable_indices;
368 std::vector<int> solution_indices;
370 indices, charges_set,
376 std::cout <<
"best_solution " << std::endl;
380 for (
Size i = 0; i < solution_indices.size(); ++i)
382 Size feature_index = variable_indices[solution_indices[i]].feature;
383 Size feature_scan_idx = variable_indices[solution_indices[i]].scan;
387 std::vector<Precursor> pcs;
389 p.
setMZ(features[feature_index].getMZ());
390 p.
setCharge(features[feature_index].getCharge());
393 ms2_spec.
setRT(scan->getRT() + rt_dist / 2.0);
400 std::cout <<
" MS2 spectra generated at: " << scan->getRT() <<
" x " << p.
getMZ() <<
"\n";
404 std::cout << solution_indices.size() <<
" out of " << features.size()
405 <<
" precursors are in best solution.\n";
411 std::cout <<
"scan based precursor selection" << std::endl;
416 std::vector<DoubleList> feature_elution_bounds;
417 std::vector<DoubleList> elution_profile_intensities;
418 std::vector<DoubleList> isotope_intensities;
420 bool meta_values_present =
false;
422 if (!features.empty() &&
423 features[0].metaValueExists(
"elution_profile_bounds") &&
424 features[0].metaValueExists(
"elution_profile_intensities") &&
425 features[0].metaValueExists(
"isotope_intensities"))
427 for (
Size feat = 0; feat < features.size(); ++feat)
429 feature_elution_bounds.push_back(features[feat].getMetaValue(
"elution_profile_bounds"));
430 elution_profile_intensities.push_back(features[feat].getMetaValue(
"elution_profile_intensities"));
431 isotope_intensities.push_back(features[feat].getMetaValue(
"isotope_intensities"));
433 meta_values_present =
true;
437 std::vector<std::vector<Size> > scan_features(experiment.
size());
439 for (
Size feat = 0; feat < features.size(); ++feat)
441 if (charges_set.count(features[feat].getCharge()))
443 Size lower_rt = features[feat].getConvexHull().getBoundingBox().minX();
444 Size upper_rt = features[feat].getConvexHull().getBoundingBox().maxX();
446 for (it = experiment.
RTBegin(lower_rt); it != experiment.
RTEnd(upper_rt); ++it)
448 scan_features[it - experiment.
begin()].push_back(feat);
453 bool dynamic_exclusion =
param_.
getValue(
"Exclusion:use_dynamic_exclusion") ==
"true" ?
true :
false;
455 ExclusionListType exclusion_list;
457 if (!dynamic_exclusion)
464 std::map<Size, typename OpenMS::DBoundingBox<2> > bounding_boxes_f;
466 for (Size feature_num = 0; feature_num < features.size(); ++feature_num)
468 if (charges_set.count(features[feature_num].getCharge()))
470 bounding_boxes_f.insert(std::make_pair(feature_num, features[feature_num].getConvexHull().getBoundingBox()));
471 const std::vector<ConvexHull2D> mass_traces = features[feature_num].getConvexHulls();
472 for (Size mass_trace_num = 0; mass_trace_num < mass_traces.size(); ++mass_trace_num)
477 bounding_boxes.insert(std::make_pair(std::make_pair(feature_num, mass_trace_num), tmp_bbox));
484 for (Size i = 0; i < experiment.
size(); ++i)
487 std::cout <<
"scan " << experiment[i].getRT() <<
":";
493 Size selected_peaks = 0, j = 0;
495 while (selected_peaks < max_spec && j < scan.size())
500 ExclusionListType::iterator it_low = exclusion_list.lower_bound(std::make_pair(peak_mz, peak_mz));
501 if (it_low != exclusion_list.end() && it_low->first.first <= peak_mz)
510 std::vector<Precursor> pcs;
511 std::set<std::pair<Size, Size> > selected_mt;
516 for (Size scan_feat_id = 0; scan_feat_id < scan_features[i].size(); ++scan_feat_id)
518 Size feature_num = scan_features[i][scan_feat_id];
519 if (bounding_boxes_f[feature_num].encloses(peak_rt, local_mz))
523 for (Size mass_trace_num = 0; mass_trace_num < features[feature_num].getConvexHulls().size(); ++mass_trace_num)
525 if (bounding_boxes[std::make_pair(feature_num, mass_trace_num)].encloses(
DPosition<2>(peak_rt, local_mz)))
527 DoubleReal elu_factor = 1.0, iso_factor = 1.0;
529 if (meta_values_present)
531 DoubleList xxx = elution_profile_intensities[feature_num];
532 DoubleList yyy = feature_elution_bounds[feature_num];
536 OPENMS_PRECONDITION(i - yyy[0] < xxx.size(),
"Tried to access invalid index for elution factor");
537 elu_factor = xxx[i - yyy[0]];
538 iso_factor = isotope_intensities[feature_num][mass_trace_num];
540 feature_intensity += features[feature_num].getIntensity() * iso_factor * elu_factor;
545 p.
setMZ(features[feature_num].getMZ());
546 p.
setCharge(features[feature_num].getCharge());
548 parent_feature_ids.push_back((
Int)feature_num);
557 ms2_spec.
setRT(experiment[i].getRT() + rt_dist / 2.0);
558 ms2_spec.
setMetaValue(
"parent_feature_ids", parent_feature_ids);
563 exclusion_list.insert(std::make_pair(std::make_pair(peak_mz - excl_window, peak_mz + excl_window), exclusion_specs + 1));
571 template <
typename InputPeakType>
575 std::vector<std::vector<std::pair<Size, Size> > > checked_mass_ranges;
577 checked_mass_ranges.reserve(mass_ranges.size());
578 for (
Size f = 0; f < mass_ranges.size(); ++f)
580 std::vector<std::pair<Size, Size> > checked_mass_ranges_f;
581 for (
Size s_idx = 0; s_idx < mass_ranges[f].size(); s_idx += 2)
583 Size s = mass_ranges[f][s_idx].first;
584 bool overlapping_features =
false;
588 const InputPeakType& peak_left_border = experiment[s][mass_ranges[f][s_idx].second];
589 const InputPeakType& peak_right_border = experiment[s][mass_ranges[f][s_idx + 1].second];
590 for (
Size fmr = 0; fmr < mass_ranges.size(); ++fmr)
594 for (
Size mr = 0; mr < mass_ranges[fmr].size(); mr += 2)
596 if (mass_ranges[fmr][mr].first == s)
598 const InputPeakType& tmp_peak_left = experiment[s][mass_ranges[fmr][mr].second];
599 const InputPeakType& tmp_peak_right = experiment[s][mass_ranges[fmr][mr + 1].second];
601 std::cout << tmp_peak_left.getMZ() <<
" < "
602 << peak_left_border.getMZ() - min_peak_distance <<
" && "
603 << tmp_peak_right.getMZ() <<
" < "
604 << peak_left_border.getMZ() - min_peak_distance <<
" ? "
605 << (tmp_peak_left.getMZ() < peak_left_border.getMZ() - min_peak_distance &&
606 tmp_peak_right.getMZ() < peak_left_border.getMZ() - min_peak_distance)
608 << tmp_peak_left.getMZ() <<
" > "
609 << peak_right_border.getMZ() + min_peak_distance <<
" && "
610 << tmp_peak_right.getMZ() <<
" > "
611 << peak_right_border.getMZ() + min_peak_distance <<
" ? "
612 << (tmp_peak_left.getMZ() > peak_right_border.getMZ() + min_peak_distance &&
613 tmp_peak_right.getMZ() > peak_right_border.getMZ() + min_peak_distance)
618 if (!((tmp_peak_left.getMZ() < peak_left_border.getMZ() - min_peak_distance &&
619 tmp_peak_right.getMZ() < peak_left_border.getMZ() - min_peak_distance) ||
620 (tmp_peak_left.getMZ() > peak_right_border.getMZ() + min_peak_distance &&
621 tmp_peak_right.getMZ() > peak_right_border.getMZ() + min_peak_distance)))
624 std::cout <<
"found overlapping peak" << std::endl;
626 overlapping_features =
true;
632 if (!overlapping_features)
635 std::cout <<
"feature in spec ok" << mass_ranges[f][s_idx].second <<
" in spec "
636 << mass_ranges[f][s_idx].first << std::endl;
638 checked_mass_ranges_f.insert(checked_mass_ranges_f.end(),
639 mass_ranges[f].begin() + s_idx,
640 mass_ranges[f].begin() + s_idx + 2);
643 checked_mass_ranges.push_back(checked_mass_ranges_f);
645 mass_ranges.swap(checked_mass_ranges);
648 template <
typename T>
651 for (
Size i = 0; i < exclusion_list.size(); ++i)
653 if (exclusion_list[i].second > 0)
654 --exclusion_list[i].second;
657 typename std::vector<std::pair<T, Size> >::iterator iter = exclusion_list.begin();
658 while (iter != exclusion_list.end() && iter->second != 0)
660 exclusion_list.erase(iter, exclusion_list.end());
667 it = exclusion_list.begin();
669 while (it != exclusion_list.end())
671 if ((it->second--) == 1)
673 exclusion_list.erase(it++);
684 #endif // OPENMS_ANALYSIS_ID_OFFLINEPRECURSORIONSELECTION_H
CoordinateType maxY() const
Accessor for max_ coordinate maximum.
Definition: DIntervalBase.h:258
A more convenient string class.
Definition: String.h:56
Precursor meta information.
Definition: Precursor.h:56
void setMinY(CoordinateType const c)
Mutator for max_ coordinate of the smaller point.
Definition: DIntervalBase.h:271
Size size() const
Definition: MSExperiment.h:117
Param param_
Container for current parameters.
Definition: DefaultParamHandler.h:148
SOLVER
Definition: LPWrapper.h:124
#define OPENMS_PRECONDITION(condition, message)
Precondition macro.
Definition: Macros.h:107
void setMaxY(CoordinateType const c)
Mutator for max_ coordinate of the larger point.
Definition: DIntervalBase.h:285
LPWrapper::SOLVER getLPSolver()
Definition: OfflinePrecursorIonSelection.h:106
CoordinateType getMZ() const
Non-mutable access to m/z.
Definition: Peak1D.h:108
ContainerType::const_iterator ConstIterator
Non-mutable iterator.
Definition: MSSpectrum.h:125
Iterator begin()
Definition: MSExperiment.h:147
bool enclosesBoundingBox(const Feature &f, typename MSExperiment< InputPeakType >::CoordinateType rt, typename MSExperiment< InputPeakType >::CoordinateType mz)
Definition: OfflinePrecursorIonSelection.h:138
void setMZ(CoordinateType mz)
Mutable access to m/z.
Definition: Peak1D.h:114
void setIntensity(IntensityType intensity)
Mutable access to the data point intensity (height)
Definition: Peak1D.h:105
ConstIterator RTEnd(CoordinateType rt) const
Fast search for spectrum range end (returns the past-the-end iterator)
Definition: MSExperiment.h:363
Class for comparison of std::pair using second ONLY e.g. for use with std::sort.
Definition: ComparatorUtils.h:368
Iterator end()
Definition: MSExperiment.h:157
ConstIterator RTBegin(CoordinateType rt) const
Fast search for spectrum range begin.
Definition: MSExperiment.h:349
CoordinateType minY() const
Accessor for max_ coordinate minimum.
Definition: DIntervalBase.h:246
const DataValue & getValue(const String &key) const
Returns a value of a parameter.
void updateExclusionList_(std::vector< std::pair< T, Size > > &exclusion_list)
Definition: OfflinePrecursorIonSelection.h:649
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:314
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:342
LPWrapper::SOLVER solver_
Definition: OfflinePrecursorIonSelection.h:134
void setMSLevel(UInt ms_level)
Sets the MS level.
Definition: MSSpectrum.h:237
An LC-MS feature.
Definition: Feature.h:66
void setLPSolver(LPWrapper::SOLVER solver)
Definition: OfflinePrecursorIonSelection.h:100
void addSpectrum(const MSSpectrum< PeakT > &spectrum)
adds a spectra to the list
Definition: MSExperiment.h:738
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:572
Representation of a mass spectrometry experiment.
Definition: MSExperiment.h:68
std::vector< SpectrumType >::const_iterator ConstIterator
Non-mutable iterator.
Definition: MSExperiment.h:103
Implements different algorithms for precursor ion selection.
Definition: OfflinePrecursorIonSelection.h:61
void setPrecursors(const std::vector< Precursor > &precursors)
sets the precursors
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:154
void setRT(DoubleReal rt)
Sets the absolute retention time (is seconds)
Definition: MSSpectrum.h:221
size_t Size
Size type e.g. used as variable which can hold result of size()
Definition: Types.h:144
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
static IntList create(const String &list)
Returns a list that is created by splitting the given comma-separated string (String are not trimmed!...
void setCharge(Int charge)
Mutable access to the charge.
A base class for all classes handling default parameters.
Definition: DefaultParamHandler.h:90
A D-dimensional bounding box.
Definition: DBoundingBox.h:51
PSLPFormulation::IndexTriple IndexTriple
Definition: OfflinePrecursorIonSelection.h:65
DoubleReal getRT() const
Definition: MSSpectrum.h:215
int Int
Signed integer type.
Definition: Types.h:100
DoubleReal list.
Definition: DoubleList.h:56
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, DoubleReal > > > &xics)
Calculate the sum of intensities of relevant features for each scan separately.
Definition: OfflinePrecursorIonSelection.h:305
double DoubleReal
Double-precision real type.
Definition: Types.h:118
Int list.
Definition: IntList.h:56