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;
119 template <
typename InputPeakType>
120 void calculateXICs_(
const FeatureMap& features,
121 const std::vector<std::vector<std::pair<Size, Size> > >& mass_ranges,
123 const std::set<Int>& charges_set,
124 std::vector<std::vector<std::pair<Size, double> > >& xics);
129 template <
typename InputPeakType>
130 void checkMassRanges_(std::vector<std::vector<std::pair<Size, Size> > >& mass_ranges,
134 void updateExclusionList_(ExclusionListType_& exclusion_list)
const;
139 template <
typename InputPeakType>
142 bool enclose_hit =
false;
144 for (
Size i = 0; i < hulls.size(); ++i)
146 if (hulls[i].getBoundingBox().encloses(rt, mz))
155 template <
typename InputPeakType>
158 std::vector<std::vector<std::pair<Size, Size> > >& indices)
160 if (experiment.
empty())
162 for (
Size f = 0; f < features.size(); ++f)
164 std::vector<std::pair<Size, Size> > vec;
166 for (
Size rt = 0; rt < experiment.
size(); ++rt)
169 if (!enclosesBoundingBox<InputPeakType>(features[f], experiment[rt].getRT(), features[f].getMZ()))
172 std::pair<Size, Size> start;
173 std::pair<Size, Size> end;
174 bool start_found =
false;
175 bool end_found =
false;
178 if (mz_iter == experiment[rt].end())
181 while (enclosesBoundingBox<InputPeakType>(features[f], experiment[rt].getRT(), mz_iter->getMZ()))
185 start.second = distance(experiment[rt].begin(), mz_iter);
186 if (mz_iter == experiment[rt].begin())
191 end.second = start.second;
194 while (mz_end != experiment[rt].end() && enclosesBoundingBox<InputPeakType>(features[f], experiment[rt].getRT(), mz_end->getMZ()))
198 end.second = distance(experiment[rt].begin(), mz_end);
201 if (start_found && end_found)
203 vec.push_back(start);
207 else if (start_found || end_found)
209 std::cout <<
"start " << start_found <<
" end " << end_found << std::endl;
210 std::cout <<
"feature: " << f <<
" rt: " << rt << std::endl;
217 std::cout << vec.size() <<
" / 2 scans" << std::endl;
218 for (
Size i = 0; i < vec.size(); i += 2)
220 std::cout <<
"Feature " << f <<
" RT : " << vec[i].first
221 <<
" MZ : " << experiment[vec[i].first][vec[i].second].getMZ() <<
" "
222 << experiment[vec[i + 1].first][vec[i + 1].second].getMZ() << std::endl;
229 std::cout <<
"According to the convex hulls no mass traces found for this feature->estimate!"
230 << features[f].getRT() <<
" " << features[f].getMZ() <<
" " << features[f].getCharge() << std::endl;
234 if (spec_iter == experiment.
end())
237 double dist1 = fabs(spec_iter->getRT() - features[f].getRT());
238 double dist2 = std::numeric_limits<double>::max();
239 double dist3 = std::numeric_limits<double>::max();
240 if ((spec_iter + 1) != experiment.
end())
242 dist2 = fabs((spec_iter + 1)->getRT() - features[f].getRT());
244 if (spec_iter != experiment.
begin())
246 dist3 = fabs((spec_iter - 1)->getRT() - features[f].getRT());
248 if (dist3 <= dist1 && dist3 <= dist2)
252 else if (dist2 <= dist3 && dist2 <= dist1)
256 std::pair<Size, Size> start;
257 std::pair<Size, Size> end;
258 start.first = distance(experiment.
begin(), spec_iter);
259 end.first = start.first;
262 if (spec_iter->
begin() == spec_iter->
end())
264 indices.push_back(vec);
267 if (mz_iter == spec_iter->
end() || (mz_iter->getMZ() > features[f].getMZ() && mz_iter != spec_iter->
begin()))
269 while (mz_iter != spec_iter->
begin())
271 if (fabs((mz_iter - 1)->getMZ() - features[f].getMZ()) < 0.5)
276 start.second = distance(spec_iter->
begin(), mz_iter);
279 std::cout << features[f].getMZ() <<
" Start: " << experiment[start.first].getRT() <<
" " << experiment[start.first][start.second].getMZ();
281 Int charge = features[f].getCharge();
284 while (mz_end + 1 != spec_iter->
end())
286 if (fabs((mz_end + 1)->getMZ() - features[f].getMZ()) < 3.0 / (
double)charge)
291 end.second = distance(spec_iter->
begin(), mz_end);
293 std::cout <<
"\tEnd: " << experiment[end.first].getRT() <<
" " << experiment[end.first][end.second].getMZ() << std::endl;
295 vec.push_back(start);
299 indices.push_back(vec);
306 template <
typename InputPeakType>
308 const std::vector<std::vector<std::pair<Size, Size> > >& mass_ranges,
310 const std::set<Int>& charges_set,
311 std::vector<std::vector<std::pair<Size, double> > >& xics)
314 xics.resize(experiment.
size());
316 for (
Size f = 0; f < mass_ranges.size(); ++f)
319 if (charges_set.count(features[f].getCharge()) < 1)
324 for (
Size s = 0; s < mass_ranges[f].size(); s += 2)
328 for (
Size j = mass_ranges[f][s].second; j <= mass_ranges[f][s + 1].second; ++j)
330 weight += experiment[mass_ranges[f][s].first][j].getIntensity();
333 xics[mass_ranges[f][s].first].push_back(std::make_pair(f, weight));
337 for (
Size s = 0; s < xics.size(); ++s)
343 template <
typename InputPeakType>
347 std::set<Int>& charges_set,
351 const double mz_isolation_window =
param_.
getValue(
"mz_isolation_window");
352 const double min_mz_peak_distance =
param_.
getValue(
"min_mz_peak_distance");
355 if (experiment.
size() > 1)
357 rt_dist = experiment[1].getRT() - experiment[0].getRT();
367 std::vector<std::vector<std::pair<Size, Size> > > indices;
370 std::vector<IndexTriple> variable_indices;
371 std::vector<int> solution_indices;
373 indices, charges_set,
379 std::cout <<
"best_solution " << std::endl;
383 for (
Size i = 0; i < solution_indices.size(); ++i)
385 Size feature_index = variable_indices[solution_indices[i]].feature;
386 Size feature_scan_idx = variable_indices[solution_indices[i]].scan;
390 std::vector<Precursor> pcs;
392 p.
setMZ(features[feature_index].getMZ());
393 p.
setCharge(features[feature_index].getCharge());
396 ms2_spec.
setRT(scan->getRT() + rt_dist / 2.0);
401 ms2_spec.
setMetaValue(
"parent_feature_ids", ListUtils::create<Int>(
String(feature_index)));
403 std::cout <<
" MS2 spectra generated at: " << scan->getRT() <<
" x " << p.
getMZ() <<
"\n";
407 std::cout << solution_indices.size() <<
" out of " << features.size()
408 <<
" precursors are in best solution.\n";
414 std::cout <<
"scan based precursor selection" << std::endl;
419 std::vector<DoubleList> feature_elution_bounds;
420 std::vector<DoubleList> elution_profile_intensities;
421 std::vector<DoubleList> isotope_intensities;
423 bool meta_values_present =
false;
425 if (!features.empty() &&
426 features[0].metaValueExists(
"elution_profile_bounds") &&
427 features[0].metaValueExists(
"elution_profile_intensities") &&
428 features[0].metaValueExists(
"isotope_intensities"))
430 for (
Size feat = 0; feat < features.size(); ++feat)
432 feature_elution_bounds.push_back(features[feat].getMetaValue(
"elution_profile_bounds"));
433 elution_profile_intensities.push_back(features[feat].getMetaValue(
"elution_profile_intensities"));
434 isotope_intensities.push_back(features[feat].getMetaValue(
"isotope_intensities"));
436 meta_values_present =
true;
440 std::vector<std::vector<Size> > scan_features(experiment.
size());
442 for (
Size feat_idx = 0; feat_idx < features.size(); ++feat_idx)
444 if (charges_set.count(features[feat_idx].getCharge()))
446 Size lower_rt = features[feat_idx].getConvexHull().getBoundingBox().minX();
447 Size upper_rt = features[feat_idx].getConvexHull().getBoundingBox().maxX();
449 for (it = experiment.
RTBegin(lower_rt); it != experiment.
RTEnd(upper_rt); ++it)
451 scan_features[it - experiment.
begin()].push_back(feat_idx);
458 int exclusion_scan_count = (int)(floor((
double)
param_.
getValue(
"Exclusion:exclusion_time") / (
double) rt_dist));
459 if (!dynamic_exclusion)
462 exclusion_scan_count = 0;
466 std::map<Size, typename OpenMS::DBoundingBox<2> > bounding_boxes_f;
468 for (
Size feature_num = 0; feature_num < features.size(); ++feature_num)
470 if (charges_set.count(features[feature_num].getCharge()))
472 bounding_boxes_f.insert(std::make_pair(feature_num, features[feature_num].getConvexHull().getBoundingBox()));
473 const std::vector<ConvexHull2D> mass_traces = features[feature_num].getConvexHulls();
474 for (
Size mass_trace_num = 0; mass_trace_num < mass_traces.size(); ++mass_trace_num)
477 tmp_bbox.
setMinY(tmp_bbox.
minY() - mz_isolation_window);
478 tmp_bbox.
setMaxY(tmp_bbox.
maxY() + mz_isolation_window);
479 bounding_boxes.insert(std::make_pair(std::make_pair(feature_num, mass_trace_num), tmp_bbox));
486 for (
Size i = 0; i < experiment.
size(); ++i)
489 std::cout <<
"scan " << experiment[i].getRT() <<
":";
496 Size selected_peaks = 0;
499 double peak_rt = scan.
getRT();
500 for (
Size j = 0; j < scan.size(); ++j)
502 if (selected_peaks >= max_spec)
break;
504 double peak_mz = scan[j].getMZ();
507 ExclusionListType_::iterator it_low = exclusion_list.lower_bound(std::make_pair(peak_mz, peak_mz));
508 if (it_low != exclusion_list.end() && it_low->first.first <= peak_mz)
517 std::vector<Precursor> pcs;
518 std::set<std::pair<Size, Size> > selected_mt;
521 double local_mz = peak_mz;
523 for (
Size scan_feat_id = 0; scan_feat_id < scan_features[i].size(); ++scan_feat_id)
525 Size feature_num = scan_features[i][scan_feat_id];
526 if (bounding_boxes_f[feature_num].encloses(peak_rt, local_mz))
529 double feature_intensity = 0;
530 for (
Size mass_trace_num = 0; mass_trace_num < features[feature_num].getConvexHulls().size(); ++mass_trace_num)
532 if (bounding_boxes[std::make_pair(feature_num, mass_trace_num)].encloses(
DPosition<2>(peak_rt, local_mz)))
534 double elu_factor = 1.0, iso_factor = 1.0;
536 if (meta_values_present)
538 DoubleList xxx = elution_profile_intensities[feature_num];
539 DoubleList yyy = feature_elution_bounds[feature_num];
543 OPENMS_PRECONDITION(i - yyy[0] < xxx.size(),
"Tried to access invalid index for elution factor");
544 elu_factor = xxx[i - yyy[0]];
545 iso_factor = isotope_intensities[feature_num][mass_trace_num];
547 feature_intensity += features[feature_num].getIntensity() * iso_factor * elu_factor;
552 p.
setMZ(features[feature_num].getMZ());
553 p.
setCharge(features[feature_num].getCharge());
555 parent_feature_ids.push_back((
Int)feature_num);
564 ms2_spec.
setRT(experiment[i].getRT() + rt_dist / 2.0);
565 ms2_spec.
setMetaValue(
"parent_feature_ids", parent_feature_ids);
570 exclusion_list.insert(std::make_pair(std::make_pair(peak_mz - min_mz_peak_distance, peak_mz + min_mz_peak_distance), exclusion_scan_count + 1));
577 template <
typename InputPeakType>
581 std::vector<std::vector<std::pair<Size, Size> > > checked_mass_ranges;
582 double min_mz_peak_distance =
param_.
getValue(
"min_mz_peak_distance");
583 checked_mass_ranges.reserve(mass_ranges.size());
584 for (
Size f = 0; f < mass_ranges.size(); ++f)
586 std::vector<std::pair<Size, Size> > checked_mass_ranges_f;
587 for (
Size s_idx = 0; s_idx < mass_ranges[f].size(); s_idx += 2)
589 Size s = mass_ranges[f][s_idx].first;
590 bool overlapping_features =
false;
594 const InputPeakType& peak_left_border = experiment[s][mass_ranges[f][s_idx].second];
595 const InputPeakType& peak_right_border = experiment[s][mass_ranges[f][s_idx + 1].second];
596 for (
Size fmr = 0; fmr < mass_ranges.size(); ++fmr)
600 for (
Size mr = 0; mr < mass_ranges[fmr].size(); mr += 2)
602 if (mass_ranges[fmr][mr].first == s)
604 const InputPeakType& tmp_peak_left = experiment[s][mass_ranges[fmr][mr].second];
605 const InputPeakType& tmp_peak_right = experiment[s][mass_ranges[fmr][mr + 1].second];
607 std::cout << tmp_peak_left.getMZ() <<
" < "
608 << peak_left_border.getMZ() - min_mz_peak_distance <<
" && "
609 << tmp_peak_right.getMZ() <<
" < "
610 << peak_left_border.getMZ() - min_mz_peak_distance <<
" ? "
611 << (tmp_peak_left.getMZ() < peak_left_border.getMZ() - min_mz_peak_distance &&
612 tmp_peak_right.getMZ() < peak_left_border.getMZ() - min_mz_peak_distance)
614 << tmp_peak_left.getMZ() <<
" > "
615 << peak_right_border.getMZ() + min_mz_peak_distance <<
" && "
616 << tmp_peak_right.getMZ() <<
" > "
617 << peak_right_border.getMZ() + min_mz_peak_distance <<
" ? "
618 << (tmp_peak_left.getMZ() > peak_right_border.getMZ() + min_mz_peak_distance &&
619 tmp_peak_right.getMZ() > peak_right_border.getMZ() + min_mz_peak_distance)
624 if (!((tmp_peak_left.getMZ() < peak_left_border.getMZ() - min_mz_peak_distance &&
625 tmp_peak_right.getMZ() < peak_left_border.getMZ() - min_mz_peak_distance) ||
626 (tmp_peak_left.getMZ() > peak_right_border.getMZ() + min_mz_peak_distance &&
627 tmp_peak_right.getMZ() > peak_right_border.getMZ() + min_mz_peak_distance)))
630 std::cout <<
"found overlapping peak" << std::endl;
632 overlapping_features =
true;
638 if (!overlapping_features)
641 std::cout <<
"feature in spec ok" << mass_ranges[f][s_idx].second <<
" in spec "
642 << mass_ranges[f][s_idx].first << std::endl;
644 checked_mass_ranges_f.insert(checked_mass_ranges_f.end(),
645 mass_ranges[f].begin() + s_idx,
646 mass_ranges[f].begin() + s_idx + 2);
649 checked_mass_ranges.push_back(checked_mass_ranges_f);
651 mass_ranges.swap(checked_mass_ranges);
657 #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
void updateExclusionList_(ExclusionListType_ &exclusion_list) const
reduce scan count for each entry, and remove every entry which has reached 0 counts ...
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:140
std::vector< Int > IntList
Vector of signed integers.
Definition: ListUtils.h:59
std::map< std::pair< double, double >, int, PairComparatorSecondElement< std::pair< double, double > > > ExclusionListType_
Definition: OfflinePrecursorIonSelection.h:114
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
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:344
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:156
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:307
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:136
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:578
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
Class for comparison of std::pair using second ONLY e.g. for use with std::sort.
Definition: ComparatorUtils.h:340
bool toBool() const
Conversion to bool.
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