35 #ifndef OPENMS_TRANSFORMATIONS_FEATUREFINDER_FEATUREFINDERALGORITHMPICKED_H
36 #define OPENMS_TRANSFORMATIONS_FEATUREFINDER_FEATUREFINDERALGORITHMPICKED_H
55 #include <boost/math/special_functions/fpclassify.hpp>
61 #include <QtCore/QDir>
83 template <
class PeakType,
class FeatureType>
117 defaults_.
setValue(
"debug",
"false",
"When debug mode is activated, several files with intermediate results are written to the folder 'debug' (do not use in parallel mode).");
120 defaults_.
setValue(
"intensity:bins", 10,
"Number of bins per dimension (RT and m/z). The higher this value, the more local the intensity significance score is.\nThis parameter should be decreased, if the algorithm is used on small regions of a map.");
122 defaults_.
setSectionDescription(
"intensity",
"Settings for the calculation of a score indicating if a peak's intensity is significant in the local environment (between 0 and 1)");
124 defaults_.
setValue(
"mass_trace:mz_tolerance", 0.03,
"Tolerated m/z deviation of peaks belonging to the same mass trace.\nIt should be larger than the m/z resolution of the instument.\nThis value must be smaller than that 1/charge_high!");
126 defaults_.
setValue(
"mass_trace:min_spectra", 10,
"Number of spectra that have to show a similar peak mass in a mass trace.");
128 defaults_.
setValue(
"mass_trace:max_missing", 1,
"Number of consecutive spectra where a high mass deviation or missing peak is acceptable.\nThis parameter should be well below 'min_spectra'!");
130 defaults_.
setValue(
"mass_trace:slope_bound", 0.1,
"The maximum slope of mass trace intensities when extending from the highest peak.\nThis parameter is important to seperate overlapping elution peaks.\nIt should be increased if feature elution profiles fluctuate a lot.");
132 defaults_.
setSectionDescription(
"mass_trace",
"Settings for the calculation of a score indicating if a peak is part of a mass trace (between 0 and 1).");
134 defaults_.
setValue(
"isotopic_pattern:charge_low", 1,
"Lowest charge to search for.");
136 defaults_.
setValue(
"isotopic_pattern:charge_high", 4,
"Highest charge to search for.");
138 defaults_.
setValue(
"isotopic_pattern:mz_tolerance", 0.03,
"Tolerated m/z deviation from the theoretical isotopic pattern.\nIt should be larger than the m/z resolution of the instument.\nThis value must be smaller than that 1/charge_high!");
140 defaults_.
setValue(
"isotopic_pattern:intensity_percentage", 10.0,
"Isotopic peaks that contribute more than this percentage to the overall isotope pattern intensity must be present.",
StringList::create(
"advanced"));
143 defaults_.
setValue(
"isotopic_pattern:intensity_percentage_optional", 0.1,
"Isotopic peaks that contribute more than this percentage to the overall isotope pattern intensity can be missing.",
StringList::create(
"advanced"));
146 defaults_.
setValue(
"isotopic_pattern:optional_fit_improvement", 2.0,
"Minimal percental improvement of isotope fit to allow leaving out an optional peak.",
StringList::create(
"advanced"));
149 defaults_.
setValue(
"isotopic_pattern:mass_window_width", 25.0,
"Window width in Dalton for precalculation of estimated isotope distributions.",
StringList::create(
"advanced"));
159 defaults_.
setSectionDescription(
"isotopic_pattern",
"Settings for the calculation of a score indicating if a peak is part of a isotopic pattern (between 0 and 1).");
161 defaults_.
setValue(
"seed:min_score", 0.8,
"Minimum seed score a peak has to reach to be used as seed.\nThe seed score is the geometric mean of intensity score, mass trace score and isotope pattern score.\nIf your features show a large deviation from the averagene isotope distribution or from an gaussian elution profile, lower this score.");
174 defaults_.
setValue(
"feature:min_score", 0.7,
"Feature score threshold for a feature to be reported.\nThe feature score is the geometric mean of the average relative deviation and the correlation between the model and the observed peaks.");
180 defaults_.
setValue(
"feature:min_trace_score", 0.5,
"Trace score threshold.\nTraces below this threshold are removed after the model fitting.\nThis parameter is important for features that overlap in m/z dimension.",
StringList::create(
"advanced"));
188 defaults_.
setValue(
"feature:rt_shape",
"symmetric",
"Choose model used for RT profile fitting. If set to symmetric a gauss shape is used, in case of asymmetric an EGH shape is used.",
StringList::create(
"advanced"));
193 defaults_.
setValue(
"feature:reported_mz",
"monoisotopic",
"The mass type that is reported for features.\n'maximum' returns the m/z value of the highest mass trace.\n'average' returns the intensity-weighted average m/z value of all contained peaks.\n'monoisotopic' returns the monoisotopic m/z value derived from the fitted isotope model.");
197 defaults_.
setValue(
"user-seed:rt_tolerance", 5.0,
"Allowed RT deviation of seeds from the user-specified seed position.");
199 defaults_.
setValue(
"user-seed:mz_tolerance", 1.1,
"Allowed m/z deviation of seeds from the user-specified seed position.");
201 defaults_.
setValue(
"user-seed:min_score", 0.5,
"Overwrites 'seed:min_score' for user-specified seeds. The cutoff is typically a bit lower in this case.");
234 Size max_isotopes = 20;
245 max_isotopes += 1000;
247 std::vector<std::pair<Size, double> > container;
248 container.push_back(std::make_pair(12, abundance_12C / 100.0));
249 container.push_back(std::make_pair(13, 1.0 - (abundance_12C / 100.0)));
250 isotopes.
set(container);
259 max_isotopes += 1000;
261 std::vector<std::pair<Size, double> > container;
262 container.push_back(std::make_pair(14, abundance_14N / 100.0));
263 container.push_back(std::make_pair(15, 1.0 - (abundance_14N / 100.0)));
264 isotopes.
set(container);
270 Param trace_fitter_params;
271 trace_fitter_params.
setValue(
"max_iteration", max_iterations);
272 trace_fitter_params.
setValue(
"epsilon_abs", epsilon_abs);
273 trace_fitter_params.
setValue(
"epsilon_rel", epsilon_rel);
279 bool user_seeds = (
seeds_.size() > 0);
289 UInt charge_count = charge_high - charge_low + 1;
293 map_[s].getFloatDataArrays().
resize(3 + 2 * charge_count);
294 map_[s].getFloatDataArrays()[0].setName(
"trace_score");
295 map_[s].getFloatDataArrays()[0].assign(scan_size, 0.0);
296 map_[s].getFloatDataArrays()[1].setName(
"intensity_score");
297 map_[s].getFloatDataArrays()[1].assign(scan_size, 0.0);
298 map_[s].getFloatDataArrays()[2].setName(
"local_max");
299 map_[s].getFloatDataArrays()[2].assign(scan_size, 0.0);
301 UInt charge = charge_low;
302 for (
Size i = 3; i < 3 + charge_count; ++i)
304 map_[s].getFloatDataArrays()[i].setName(
String(
"pattern_score_") + charge);
305 map_[s].getFloatDataArrays()[i].assign(scan_size, 0.0);
310 for (
Size i = 3 + charge_count; i < 3 + 2 * charge_count; ++i)
312 map_[s].getFloatDataArrays()[i].setName(
String(
"overall_score_") + charge);
313 map_[s].getFloatDataArrays()[i].assign(scan_size, 0.0);
324 dir.mkpath(
"debug/features");
325 log_.open(
"debug/log.txt");
332 if (
debug_)
log_ <<
"Precalculating intensity thresholds ..." << std::endl;
345 DoubleReal max_rt = rt_start + (rt + 1) * intensity_rt_step_;
346 std::vector<DoubleReal> tmp;
351 DoubleReal max_mz = mz_start + (mz + 1) * intensity_mz_step_;
357 tmp.push_back(it->getIntensity());
364 std::sort(tmp.begin(), tmp.end());
365 for (
Size i = 0; i < 21; ++i)
367 Size index = (
Size) std::floor(0.05 * i * (tmp.size() - 1));
392 ff_->
startProgress(min_spectra_, end_iteration,
"Precalculating mass trace scores");
394 for (
Size s = min_spectra_; s < end_iteration; ++s)
399 for (
Size p = 0; p < spectrum.size(); ++p)
401 std::vector<DoubleReal> scores;
402 scores.reserve(2 * min_spectra_);
405 Real inte = spectrum[p].getIntensity();
408 bool is_max_peak =
true;
413 Size spec_index =
map_[s + i].findNearest(pos);
415 if (position_score > 0 &&
map_[s + i][spec_index].getIntensity() > inte) is_max_peak =
false;
416 scores.push_back(position_score);
420 scores.push_back(0.0);
427 Size spec_index =
map_[s - i].findNearest(pos);
429 if (position_score > 0 &&
map_[s - i][spec_index].getIntensity() > inte) is_max_peak =
false;
430 scores.push_back(position_score);
434 scores.push_back(0.0);
438 DoubleReal trace_score = std::accumulate(scores.begin(), scores.end(), 0.0) / scores.size();
441 map_[s].getFloatDataArrays()[0][p] = trace_score;
442 map_[s].getFloatDataArrays()[2][p] = is_max_peak;
456 ff_->
startProgress(0, num_isotopes,
"Precalculating isotope distributions");
462 for (
Size index = 0; index < num_isotopes; ++index)
483 bool is_begin =
true;
490 if (!is_end && !is_begin) is_end =
true;
491 if (is_begin) ++begin;
492 else if (is_end) ++end;
526 Int plot_nr_global = -1;
527 Int feature_nr_global = 0;
530 UInt meta_index_isotope = 3 +
c - charge_low;
531 UInt meta_index_overall = 3 + charge_count +
c - charge_low;
533 Size feature_candidates = 0;
534 std::vector<Seed> seeds;
544 for (
Size p = 0; p < spectrum.size(); ++p)
556 for (
Size i = 0; i < isotopes.
size(); ++i)
565 if (pattern_score > 0.0)
567 for (
Size i = 0; i < pattern.peak.size(); ++i)
569 if (pattern.peak[i] >= 0 && pattern_score >
map_[pattern.spectrum[i]].getFloatDataArrays()[meta_index_isotope][pattern.peak[i]])
571 map_[pattern.spectrum[i]].getFloatDataArrays()[meta_index_isotope][pattern.peak[i]] = pattern_score;
587 for (
Size s = min_spectra_; s < end_of_iteration; ++s)
595 DoubleReal overall_score = std::pow(meta[0][p] * meta[1][p] * meta[meta_index_isotope][p], 1.0f / 3.0f);
596 meta[meta_index_overall][p] = overall_score;
599 if (meta[2][p] != 0.0)
602 if (!user_seeds && overall_score >= min_seed_score)
608 seeds.push_back(seed);
611 else if (user_seeds && overall_score >= user_seed_score)
615 tmp.setMZ(
map_[s][p].getMZ() - user_mz_tol);
616 for (
typename FeatureMapType::const_iterator it = std::lower_bound(
seeds_.begin(),
seeds_.end(), tmp,
typename FeatureType::MZLess()); it <
seeds_.end(); ++it)
618 if (it->getMZ() >
map_[s][p].getMZ() + user_mz_tol)
622 if (fabs(it->getMZ() -
map_[s][p].getMZ()) < user_mz_tol &&
623 fabs(it->getRT() -
map_[s].getRT()) < user_rt_tol)
629 seeds.push_back(seed);
638 std::sort(seeds.rbegin(), seeds.rend());
644 seed_map.reserve(seeds.size());
645 for (
Size i = 0; i < seeds.size(); ++i)
647 Size spectrum = seeds[i].spectrum;
648 Size peak = seeds[i].peak;
656 tmp.
setMetaValue(
"pattern_score", meta[meta_index_isotope][peak]);
658 seed_map.push_back(tmp);
664 std::cout <<
"Found " << seeds.size() <<
" seeds for charge " <<
c <<
"." << std::endl;
679 std::map<Size, std::vector<Size> > seeds_in_features;
681 FeatureMapType tmp_feature_map;
685 #pragma omp parallel for
695 const PeakType& peak = spectrum[seeds[i].peak];
703 log_ << std::endl <<
"Seed " << i <<
":" << std::endl;
706 log_ <<
" - RT: " << spectrum.
getRT() << std::endl;
707 log_ <<
" - MZ: " << peak.
getMZ() << std::endl;
718 abort_(seeds[i],
"Could not find good enough isotope pattern containing the seed");
726 traces.reserve(best_pattern.
peak.size());
730 DoubleReal seed_mz =
map_[seeds[i].spectrum][seeds[i].peak].getMZ();
734 abort_(seeds[i],
"Could not extend seed");
747 #pragma omp critical (FeatureFinderAlgorithmPicked_PLOTNR)
750 plot_nr = ++plot_nr_global;
760 traces[traces.
max_trace].updateMaximum();
763 double egh_tau = 0.0;
772 alt_p.setValue(
"max_iteration", max_iterations);
773 alt_p.setValue(
"epsilon_abs", epsilon_abs);
774 alt_p.setValue(
"epsilon_rel", epsilon_rel);
777 alt_fitter->
fit(traces);
807 bool feature_ok =
checkFeatureQuality_(fitter, new_traces, seed_mz, min_feature_score, error_msg, fit_score, correlation, final_score);
809 #pragma omp critical (FeatureFinderAlgorithmPicked_DEBUG)
824 abort_(seeds[i], error_msg);
862 for (
Size t = 0; t < traces.size(); ++t)
864 for (
Size p = 0; p < traces[t].peaks.size(); ++p)
866 average_mz += traces[t].peaks[p].second->getMZ() * traces[t].peaks[p].second->getIntensity();
867 total_intensity += traces[t].peaks[p].second->getIntensity();
870 average_mz /= total_intensity;
891 for (
Size j = 0; j < traces.size(); ++j)
897 #pragma omp critical (FeatureFinderAlgorithmPicked_TMPFEATUREMAP)
900 tmp_feature_map[i] = f;
906 for (
Size j = i + 1; j < seeds.size(); ++j)
913 #pragma omp critical (FeatureFinderAlgorithmPicked_SEEDSINFEATURES)
916 seeds_in_features[i].push_back(j);
929 std::vector<Size> seeds_contained;
930 for (
typename std::map<Size, FeatureType>::iterator iter = tmp_feature_map.begin(); iter != tmp_feature_map.end(); ++iter)
932 Size seed_nr = iter->first;
933 bool is_used =
false;
934 for (
Size i = 0; i < seeds_contained.size(); ++i)
936 if (seed_nr == seeds_contained[i]) { is_used =
true;
break; }
940 ++feature_candidates;
943 iter->second.setMetaValue(3, feature_nr_global);
947 std::vector<Size> curr_seed = seeds_in_features[seed_nr];
948 for (
Size k = 0;
k < curr_seed.size(); ++
k)
950 seeds_contained.push_back(curr_seed[
k]);
957 std::cout <<
"Found " << feature_candidates <<
" feature candidates for charge " <<
c <<
"." << std::endl;
966 if (
debug_)
log_ <<
"Resolving intersecting features (" <<
features_->size() <<
" candidates)" << std::endl;
970 std::vector<DBoundingBox<2> > bbs(
features_->size());
975 bbs[i] = (*features_)[i].getConvexHull().getBoundingBox();
976 if (bbs[i].height() > max_mz_span)
978 max_mz_span = bbs[i].height();
992 if (f2.
getMZ() - f1.
getMZ() > 2.0 * max_mz_span)
break;
996 if (!bbs[i].intersects(bbs[j]))
continue;
1004 if (
debug_)
log_ <<
" - Intersection (" << (i + 1) <<
"/" << (j + 1) <<
"): " << intersection << std::endl;
1009 if (
debug_)
log_ <<
" - same charge -> removing duplicate " << (j + 1) << std::endl;
1015 if (
debug_)
log_ <<
" - same charge -> removing duplicate " << (i + 1) << std::endl;
1022 if (
debug_)
log_ <<
" - different charge (one is the multiple of the other) -> removing lower charge " << (i + 1) << std::endl;
1028 if (
debug_)
log_ <<
" - different charge (one is the multiple of the other) -> removing lower charge " << (i + 1) << std::endl;
1036 if (
debug_)
log_ <<
" - different charge -> removing lower score " << (j + 1) << std::endl;
1042 if (
debug_)
log_ <<
" - different charge -> removing lower score " << (i + 1) << std::endl;
1050 LOG_INFO <<
"Removed " << removed <<
" overlapping features." << std::endl;
1056 if (
features_->operator[](i).getIntensity() != 0.0)
1058 tmp.push_back(
features_->operator[](i));
1065 std::cout <<
features_->size() <<
" features left." << std::endl;
1068 std::cout << std::endl;
1069 std::cout <<
"Abort reasons during feature construction:" << std::endl;
1070 for (std::map<String, UInt>::const_iterator it =
aborts_.begin(); it !=
aborts_.end(); ++it)
1072 std::cout <<
"- " << it->first <<
": " << it->second << std::endl;
1083 f.
setRT(
map_[it2->first.spectrum].getRT());
1084 f.
setMZ(
map_[it2->first.spectrum][it2->first.peak].getMZ());
1088 abort_map.push_back(f);
1096 map_[s].getFloatDataArrays().erase(
map_[s].getFloatDataArrays().begin() + 2);
1110 return "centroided";
1184 if (
debug_)
log_ <<
"Abort: " << reason << std::endl;
1198 for (
Size i = 0; i < hulls1.size(); ++i)
1200 s1 += hulls1[i].getBoundingBox().width();
1206 for (
Size j = 0; j < hulls2.size(); ++j)
1208 s2 += hulls2[j].getBoundingBox().width();
1213 for (
Size i = 0; i < hulls1.size(); ++i)
1216 for (
Size j = 0; j < hulls2.size(); ++j)
1224 overlap += bb2.
width();
1229 overlap += bb1.
width();
1245 return overlap / std::min(s1, s2);
1272 if (
debug_)
log_ <<
"Testing isotope patterns for charge " << charge <<
": " << std::endl;
1275 if (
debug_)
log_ <<
" - Seed: " << center.
peak <<
" (mz:" << spectrum[center.
peak].getMZ() <<
")" << std::endl;
1279 if (
debug_)
log_ <<
" - Mass window: " << mass_window << std::endl;
1281 while (end < spectrum.size() &&
1282 spectrum[end].getMZ() < spectrum[center.
peak].getMZ() + mass_window)
1290 while (begin >= 0 &&
1291 spectrum[begin].getMZ() > spectrum[center.
peak].getMZ() - mass_window)
1296 if (
debug_)
log_ <<
" - Begin: " << begin <<
" (mz:" << spectrum[begin].getMZ() <<
")" << std::endl;
1297 if (
debug_)
log_ <<
" - End: " << end <<
" (mz:" << spectrum[end].getMZ() <<
")" << std::endl;
1301 for (
Size start = begin; start <= end; ++start)
1304 Size peak_index = start;
1306 if (
debug_)
log_ <<
" - Fitting at " << start <<
" (mz:" << spectrum[start].getMZ() <<
")" << std::endl;
1307 for (
Size iso = 0; iso < isotopes.
size(); ++iso)
1314 bool seed_contained =
false;
1315 for (
Size iso = 0; iso < pattern.peak.size(); ++iso)
1317 if (pattern.peak[iso] == (
Int)center.
peak && pattern.spectrum[iso] == center.
spectrum)
1319 seed_contained =
true;
1323 if (!seed_contained)
1325 if (
debug_)
log_ <<
" - aborting: seed is not contained!" << std::endl;
1332 seed_contained =
false;
1333 for (
Size iso = 0; iso < pattern.peak.size(); ++iso)
1335 if (pattern.peak[iso] == (
Int)center.
peak &&
1336 pattern.spectrum[iso] == center.
spectrum)
1338 seed_contained =
true;
1342 if (!seed_contained)
1344 if (
debug_)
log_ <<
" - aborting: seed was removed during isotope fit!" << std::endl;
1348 if (
debug_)
log_ <<
" - final score: " << score << std::endl;
1349 if (score > max_score)
1352 best_pattern = pattern;
1355 if (
debug_)
log_ <<
" - best score : " << max_score << std::endl;
1371 Size max_trace_index = 0;
1372 for (
Size p = 0; p < pattern.
peak.size(); ++p)
1374 if (pattern.
peak[p] < 0)
continue;
1378 max_trace_index = p;
1387 if (
debug_)
log_ <<
" - Trace " << max_trace_index <<
" (maximum intensity)" << std::endl;
1388 if (
debug_)
log_ <<
" - extending from: " <<
map_[start_index].getRT() <<
" / " << start_mz <<
" (int: " << start_peak->
getIntensity() <<
")" << std::endl;
1391 max_trace.
peaks.push_back(std::make_pair(start_rt, start_peak));
1392 extendMassTrace_(max_trace, start_index, start_mz,
false, meta_index_overall);
1393 extendMassTrace_(max_trace, start_index, start_mz,
true, meta_index_overall);
1397 if (
debug_)
log_ <<
" - rt bounds: " << rt_min <<
"-" << rt_max << std::endl;
1401 if (
debug_)
log_ <<
" - could not extend trace with maximum intensity => abort" << std::endl;
1404 for (
Size p = 0; p < pattern.
peak.size(); ++p)
1406 if (
debug_)
log_ <<
" - Trace " << p << std::endl;
1407 if (p == max_trace_index)
1409 if (
debug_)
log_ <<
" - previously extended maximum trace" << std::endl;
1410 traces.push_back(max_trace);
1417 starting_peak.
peak = pattern.
peak[p];
1418 if (pattern.
peak[p] == -2)
1420 if (
debug_)
log_ <<
" - removed during isotope fit" << std::endl;
1423 else if (pattern.
peak[p] == -1)
1425 if (
debug_)
log_ <<
" - missing" << std::endl;
1436 for (
Size spectrum_index = begin; spectrum_index < end; ++spectrum_index)
1442 peak_index =
map_[spectrum_index].findNearest(
map_[starting_peak.
spectrum][starting_peak.
peak].getMZ());
1449 if (peak_index < 0 ||
1450 map_[spectrum_index][peak_index].getIntensity() <= inte ||
1457 starting_peak.
spectrum = spectrum_index;
1458 starting_peak.
peak = peak_index;
1459 inte =
map_[spectrum_index][peak_index].getIntensity();
1468 trace.
peaks.push_back(std::make_pair(
map_[starting_peak.
spectrum].getRT(), seed));
1475 if (
debug_)
log_ <<
" - could not extend trace " << std::endl;
1487 traces.push_back(trace);
1516 std::reverse(trace.
peaks.begin(), trace.
peaks.end());
1524 bool boundaries =
false;
1525 if (max_rt != min_rt)
1534 std::vector<DoubleReal> deltas(delta_count - 1, 0);
1536 DoubleReal last_observed_intensity = trace.
peaks.back().second->getIntensity();
1538 UInt missing_peaks = 0;
1539 Size peaks_before_extension = trace.
peaks.size();
1540 String abort_reason =
"";
1542 while ((!increase_rt && spectrum_index >= 0) || (increase_rt && spectrum_index < (
SignedSize)
map_.
size()))
1545 ((!increase_rt &&
map_[spectrum_index].getRT() < min_rt) ||
1546 (increase_rt &&
map_[spectrum_index].getRT() > max_rt))
1549 abort_reason =
"Hit upper/lower boundary";
1557 peak_index =
map_[spectrum_index].findNearest(mz);
1567 ||
map_[spectrum_index].getFloatDataArrays()[meta_index_overall][peak_index] < 0.01
1575 abort_reason =
"too many peaks missing";
1584 trace.
peaks.push_back(std::make_pair(
map_[spectrum_index].getRT(), &(
map_[spectrum_index][peak_index])));
1587 deltas.push_back((
map_[spectrum_index][peak_index].getIntensity() - last_observed_intensity) / last_observed_intensity);
1588 last_observed_intensity =
map_[spectrum_index][peak_index].getIntensity();
1591 DoubleReal average_delta = std::accumulate(deltas.end() - delta_count, deltas.end(), 0.0) / (
DoubleReal)delta_count;
1592 if (average_delta > current_slope_bound)
1594 abort_reason =
String(
"Average delta above threshold: ") + average_delta +
"/" + current_slope_bound;
1597 Size remove = std::min((
Size)(trace.
peaks.size() - peaks_before_extension), delta_count - 1);
1604 if (increase_rt) ++spectrum_index;
1605 else --spectrum_index;
1607 if (
debug_)
log_ <<
" - Added " << (trace.
peaks.size() - peaks_before_extension) <<
" peaks (abort: " << abort_reason <<
")" << std::endl;
1611 template <
typename SpectrumType>
1615 DoubleReal distance = std::fabs(pos - spec[index].getMZ());
1617 while (index < spec.size())
1619 DoubleReal new_distance = std::fabs(pos - spec[index].getMZ());
1620 if (new_distance < distance)
1622 distance = new_distance;
1644 if (
debug_)
log_ <<
" - Isotope " << pattern_index <<
": ";
1652 peak_index =
nearest_(pos, spectrum, peak_index);
1656 if (mz_score != 0.0)
1659 pattern.
peak[pattern_index] = peak_index;
1660 pattern.
spectrum[pattern_index] = spectrum_index;
1661 intensity += spectrum[peak_index].getIntensity();
1662 pos_score += mz_score;
1667 if (spectrum_index != 0 && !
map_[spectrum_index - 1].empty())
1672 if (mz_score != 0.0)
1675 intensity += spectrum_before[index_before].getIntensity();
1676 pos_score += mz_score;
1679 if (pattern.
peak[pattern_index] == -1)
1681 pattern.
peak[pattern_index] = index_before;
1682 pattern.
spectrum[pattern_index] = spectrum_index - 1;
1692 if (mz_score != 0.0)
1695 intensity += spectrum_after[index_after].getIntensity();
1696 pos_score += mz_score;
1699 if (pattern.
peak[pattern_index] == -1)
1701 pattern.
peak[pattern_index] = index_after;
1702 pattern.
spectrum[pattern_index] = spectrum_index + 1;
1710 pattern.
peak[pattern_index] = -1;
1711 pattern.
mz_score[pattern_index] = 0.0;
1716 if (
debug_)
log_ <<
"=> " << intensity / matches << std::endl;
1717 pattern.
mz_score[pattern_index] = pos_score / matches;
1718 pattern.
intensity[pattern_index] = intensity / matches;
1726 if (diff <= 0.5 * allowed_deviation)
1728 return 0.1 * (0.5 * allowed_deviation - diff) / (0.5 * allowed_deviation) + 0.9;
1730 else if (diff <= allowed_deviation)
1732 return 0.9 * (allowed_deviation - diff) / (0.5 * allowed_deviation);
1744 if (pattern.
peak[iso] == -1)
1746 if (
debug_)
log_ <<
" - aborting: core peak is missing" << std::endl;
1754 Size best_begin = 0;
1757 if (pattern.
peak[i - 1] == -1)
1766 if (pattern.
peak[pattern.
peak.size() - i] == -1)
1772 if (
debug_)
log_ <<
" - best_begin/end: " << best_begin <<
"/" << best_end << std::endl;
1778 if (isotopes.
size() - b - e > 2 || (b == best_begin &&
1780 isotopes.
size() - b - e > 1))
1783 if (boost::math::isnan(int_score)) int_score = 0.0;
1785 if (
debug_)
log_ <<
" - fit (" << b <<
"/" << e <<
"): " << int_score;
1789 best_int_score = int_score;
1799 if (pattern.
mz_score.size() - best_begin - best_end == 0)
1805 for (
Size i = 0; i < best_begin; ++i)
1807 pattern.
peak[i] = -2;
1812 for (
Size i = 0; i < best_end; ++i)
1814 pattern.
peak[isotopes.
size() - 1 - i] = -2;
1819 if (consider_mz_distances)
1821 best_int_score *= std::accumulate(pattern.
mz_score.begin() + best_begin, pattern.
mz_score.end() - best_end, 0.0) / (pattern.
mz_score.size() - best_begin - best_end);
1825 OPENMS_POSTCONDITION(best_int_score >= 0.0, (
String(
"Internal error: Isotope score (") + best_int_score +
") should be >=0.0").c_str())
1826 OPENMS_POSTCONDITION(best_int_score <= 1.0, (
String(
"Internal error: Isotope score (") + best_int_score +
") should be <=1.0").c_str())
1827 return best_int_score;
1860 mh = mz_bin / 2 + 1;
1864 ml = mz_bin / 2 - 1;
1877 rh = rt_bin / 2 + 1;
1881 rl = rt_bin / 2 - 1;
1891 DoubleReal d1 = std::sqrt(std::pow(1.0 - drl, 2) + std::pow(1.0 - dml, 2));
1892 DoubleReal d2 = std::sqrt(std::pow(1.0 - drh, 2) + std::pow(1.0 - dml, 2));
1893 DoubleReal d3 = std::sqrt(std::pow(1.0 - drl, 2) + std::pow(1.0 - dmh, 2));
1894 DoubleReal d4 = std::sqrt(std::pow(1.0 - drh, 2) + std::pow(1.0 - dmh, 2));
1919 LOG_DEBUG <<
"use asymmetric rt peak shape" << std::endl;
1925 LOG_DEBUG <<
"use symmetric rt peak shape" << std::endl;
1935 std::vector<DoubleReal>::const_iterator it = std::lower_bound(quantiles20.begin(), quantiles20.end(), intensity);
1937 if (it == quantiles20.end())
1943 if (it == quantiles20.begin())
1945 bin_score = 0.05 * intensity / *it;
1950 bin_score = 0.05 * (intensity - *(it - 1)) / (*it - *(it - 1));
1954 0.05 * ((it - quantiles20.begin()) - 1.0);
1957 if (
final < 0.0)
final = 0.0;
1958 if (
final > 1.0)
final = 1.0;
1986 if (
debug_)
log_ <<
" => RT bounds: " << low_bound <<
" - " << high_bound << std::endl;
1987 for (
Size t = 0; t < traces.size(); ++t)
1995 std::vector<DoubleReal> v_theo, v_real;
1999 if (trace.
peaks[
k].first >= low_bound && trace.
peaks[
k].first <= high_bound)
2005 v_theo.push_back(theo);
2007 v_real.push_back(real);
2008 deviation += std::fabs(real - theo) / theo;
2014 if (!new_trace.
peaks.empty())
2016 fit_score = deviation / new_trace.
peaks.size();
2018 final_score = std::sqrt(correlation * std::max(0.0, 1.0 - fit_score));
2020 if (
debug_)
log_ <<
" - peaks: " << new_trace.
peaks.size() <<
" / " << trace.
peaks.size() <<
" - relative deviation: " << fit_score <<
" - correlation: " << correlation <<
" - final score: " << correlation << std::endl;
2027 if (
debug_)
log_ <<
" - removed this and previous traces due to bad fit" << std::endl;
2034 if (
debug_)
log_ <<
" - aborting (max trace was removed)" << std::endl;
2039 if (
debug_)
log_ <<
" - removed due to bad fit => omitting the rest" << std::endl;
2047 new_traces.push_back(new_trace);
2050 new_traces.
max_trace = new_traces.size() - 1;
2085 bool feature_ok =
true;
2094 error_msg =
"Invalid fit: Fitted model is bigger than 'max_rt_span'";
2102 error_msg =
"Invalid feature after fit - too few traces or peaks left";
2108 std::pair<DoubleReal, DoubleReal> rt_bounds = feature_traces.
getRTBounds();
2112 error_msg =
"Invalid fit: Center outside of feature bounds";
2119 std::pair<DoubleReal, DoubleReal> rt_bounds = feature_traces.
getRTBounds();
2123 error_msg =
"Invalid fit: Less than 'min_rt_span' left after fit";
2130 std::vector<DoubleReal> v_theo, v_real;
2132 for (
Size t = 0; t < feature_traces.size(); ++t)
2139 v_theo.push_back(theo);
2141 v_real.push_back(real);
2142 deviation += std::fabs(real - theo) / theo;
2145 fit_score = std::max(0.0, 1.0 - (deviation / feature_traces.
getPeakCount()));
2147 final_score = std::sqrt(correlation * fit_score);
2149 if (final_score < min_feature_score)
2152 error_msg =
"Feature quality too low after fit";
2158 log_ <<
"Quality estimation:" << std::endl;
2159 log_ <<
" - relative deviation: " << fit_score << std::endl;
2160 log_ <<
" - correlation: " << correlation << std::endl;
2161 log_ <<
" => final score: " << final_score << std::endl;
2185 const String path =
"debug/features/")
2193 for (
Size k = 0;
k < traces.size(); ++
k)
2195 for (
Size j = 0; j < traces[
k].peaks.size(); ++j)
2197 tf.push_back(
String(pseudo_rt_shift *
k + traces[
k].peaks[j].first) +
"\t" + traces[
k].peaks[j].second->getIntensity());
2200 tf.
store(path + plot_nr +
".dta");
2205 for (
Size k = 0;
k < new_traces.size(); ++
k)
2207 for (
Size j = 0; j < new_traces[
k].peaks.size(); ++j)
2209 tf.push_back(
String(pseudo_rt_shift *
k + new_traces[
k].peaks[j].first) +
"\t" + new_traces[
k].peaks[j].second->getIntensity());
2213 tf.
store(path + plot_nr +
"_cropped.dta");
2214 script = script +
", \"" + path + plot_nr +
"_cropped.dta\" title 'feature ";
2218 script = script +
" - " + error_msg;
2224 script = script +
"' with points 3";
2228 for (
Size k = 0;
k < traces.size(); ++
k)
2234 script = script +
", " + fun +
"(x) title 'Trace " + k +
" (m/z: " +
String::number(traces[k].getAvgMZ(), 4) +
")'";
2238 tf.push_back(
"set xlabel \"pseudo RT (mass traces side-by-side)\"");
2239 tf.push_back(
"set ylabel \"intensity\"");
2240 tf.push_back(
"set samples 1000");
2241 tf.push_back(script);
2242 tf.push_back(
"pause -1");
2243 tf.
store(path + plot_nr +
".plot");
2257 #endif // OPENMS_TRANSFORMATIONS_FEATUREFINDER_FEATUREFINDERALGORITHMPICKED_H
void sortByMZ()
Sort features by m/z position.
Definition: FeatureMap.h:322
bool encloses(const PositionType &position) const
Checks whether this range contains a certain point.
Definition: DBoundingBox.h:159
std::vector< std::vector< std::vector< DoubleReal > > > intensity_thresholds_
Precalculated intensity 20-quantiles (binned)
Definition: FeatureFinderAlgorithmPicked.h:1154
bool debug_
debug flag
Definition: FeatureFinderAlgorithmPicked.h:1119
bool isOdd(UInt x)
Returns true if the given interger is odd.
Definition: MathFunctions.h:124
virtual DoubleReal getFeatureIntensityContribution()=0
MapType::SpectrumType SpectrumType
Definition: FeatureFinderAlgorithmPicked.h:93
const ChargeType & getCharge() const
Non-mutable access to charge state.
Size size() const
returns the size of the distribtion which is the number of isotopes in the distribution ...
virtual String getGnuplotFormula(FeatureFinderAlgorithmPickedHelperStructs::MassTrace< PeakType > const &trace, const char function_name, const DoubleReal baseline, const DoubleReal rt_shift)=0
ContainerType::iterator Iterator
Definition: IsotopeDistribution.h:70
FeatureMapType * features_
Output data pointer.
Definition: FeatureFinderAlgorithm.h:144
float Real
Real type.
Definition: Types.h:109
Param defaults_
Container for default parameters. This member should be filled in the constructor of derived classes!...
Definition: DefaultParamHandler.h:155
DoubleReal positionScore_(DoubleReal pos1, DoubleReal pos2, DoubleReal allowed_deviation) const
Calculates a score between 0 and 1 for the m/z deviation of two peaks.
Definition: FeatureFinderAlgorithmPicked.h:1723
DoubleReal intersection_(const Feature &f1, const Feature &f2) const
Definition: FeatureFinderAlgorithmPicked.h:1193
void setValue(const String &key, const DataValue &value, const String &description="", const StringList &tags=StringList())
Sets a value.
TraceFitter< PeakType > * chooseTraceFitter_(double &tau)
Choose a the best trace fitter for the current mass traces based on the user parameter (symmetric...
Definition: FeatureFinderAlgorithmPicked.h:1914
A more convenient string class.
Definition: String.h:56
DoubleReal intensity_percentage_
Isotope pattern intensity contribution of required peaks.
Definition: FeatureFinderAlgorithmPicked.h:1134
DoubleReal intensityScore_(Size spectrum, Size peak) const
Compute the intensity score for the peak peak in spectrum spectrum.
Definition: FeatureFinderAlgorithmPicked.h:1840
Size spectrum
Spectrum index.
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:65
DoubleReal min_isotope_fit_
Mimimum isotope pattern fit for a feature.
Definition: FeatureFinderAlgorithmPicked.h:1139
FeatureFinder * ff_
Pointer to the calling FeatureFinder that is used to access the feature flags.
Definition: FeatureFinderAlgorithm.h:147
A 2-dimensional raw data point or peak.
Definition: Peak2D.h:55
Size size() const
Definition: MSExperiment.h:117
void setMZ(CoordinateType coordinate)
Mutable access to the m/z coordinate (index 1)
Definition: Peak2D.h:197
FeatureFinderAlgorithm for picked peaks.
Definition: FeatureFinderAlgorithmPicked.h:84
QualityType getOverallQuality() const
Non-mutable access to the overall quality.
#define LOG_INFO
Macro if a information, e.g. a status should be reported.
Definition: LogStream.h:455
void cropFeature_(TraceFitter< PeakType > *fitter, const MassTraces &traces, MassTraces &new_traces)
Creates new mass traces new_traces based on the fitting result and the original traces traces...
Definition: FeatureFinderAlgorithmPicked.h:1979
Size max_trace
Maximum intensity trace.
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:269
Size findNearest(CoordinateType mz) const
Binary search for the peak nearest to a specific m/z.
Definition: MSSpectrum.h:476
void extendMassTraces_(const IsotopePattern &pattern, MassTraces &traces, Size meta_index_overall) const
Definition: FeatureFinderAlgorithmPicked.h:1367
Param param_
Container for current parameters.
Definition: DefaultParamHandler.h:148
ConvexHull2D & getConvexHull() const
Returns the overall convex hull of the feature (calculated from the convex hulls of the mass traces) ...
void findIsotope_(DoubleReal pos, Size spectrum_index, IsotopePattern &pattern, Size pattern_index, Size &peak_index) const
Searches for an isotopic peak in the current spectrum and the adjacent spectra.
Definition: FeatureFinderAlgorithmPicked.h:1642
IntensityType getIntensity() const
Definition: Peak2D.h:161
DoubleReal mass_window_width_
Width of the isotope pattern mass bins.
Definition: FeatureFinderAlgorithmPicked.h:1137
Helper struct for a collection of mass traces used in FeatureFinderAlgorithmPicked.
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:153
std::ofstream log_
Output stream for log/debug info.
Definition: FeatureFinderAlgorithmPicked.h:1117
DoubleReal pattern_tolerance_
Stores mass_trace:mz_tolerance.
Definition: FeatureFinderAlgorithmPicked.h:1129
void store(const String &filename, const FeatureMap<> &feature_map)
stores the map feature_map in file with name filename.
Abstract fitter for RT profile fitting.
Definition: TraceFitter.h:62
DoubleReal min_rt_span_
Minimum RT range that has to be left after the fit.
Definition: FeatureFinderAlgorithmPicked.h:1141
bool isValid(DoubleReal seed_mz, DoubleReal trace_tolerance)
Checks if still valid (seed still contained and enough traces)
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:174
A container for features.
Definition: FeatureMap.h:111
Iterator end()
Definition: IsotopeDistribution.h:187
#define OPENMS_POSTCONDITION(condition, message)
Postcondition macro.
Definition: Macros.h:114
unsigned int UInt
Unsigned integer type.
Definition: Types.h:92
virtual bool checkMinimalRTSpan(const std::pair< DoubleReal, DoubleReal > &rt_bounds, const DoubleReal min_rt_span)=0
Size optional_begin
Number of optional peaks at the beginning of the pattern.
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:282
Helper structure for a theoretical isotope pattern used in FeatureFinderAlgorithmPicked.
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:277
Abstract base class for FeatureFinder algorithms.
Definition: FeatureFinderAlgorithm.h:74
const TheoreticalIsotopePattern & getIsotopeDistribution_(DoubleReal mass) const
Returns the isotope distribution for a certain mass window.
Definition: FeatureFinderAlgorithmPicked.h:1249
std::vector< DoubleReal > mz_score
m/z score of peak (0 if peak index is -1 or -2)
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:309
Forward iterator for an area of peaks in an experiment.
Definition: AreaIterator.h:58
Isotope distribution class.
Definition: IsotopeDistribution.h:61
Size size() const
Returns the size.
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:290
std::map< Seed, String > abort_reasons_
Array of abort reasons.
Definition: FeatureFinderAlgorithmPicked.h:1123
FeatureFinderAlgorithm< PeakType, FeatureType >::FeatureMapType FeatureMapType
Definition: FeatureFinderAlgorithmPicked.h:92
bool encloses(DoubleReal rt, DoubleReal mz) const
Returns if the mass trace convex hulls of the feature enclose the position specified by rt and mz...
Size setUniqueId()
Assigns a new, valid unique id. Always returns 1.
Definition: UniqueIdInterface.h:150
static String number(DoubleReal d, UInt n)
returns a string for d with exactly n decimal places
virtual DoubleReal getLowerRTBound() const =0
void setSectionDescription(const String &key, const String &description)
Sets a description for an existing section.
String reported_mz_
The mass type that is reported for features. 'maximum' returns the m/z value of the highest mass trac...
Definition: FeatureFinderAlgorithmPicked.h:1144
DoubleReal trace_tolerance_
Stores isotopic_pattern:mz_tolerance.
Definition: FeatureFinderAlgorithmPicked.h:1130
ptrdiff_t SignedSize
Signed Size type e.g. used as pointer difference.
Definition: Types.h:151
bool intersects(const DBoundingBox &bounding_box) const
Definition: DBoundingBox.h:180
std::vector< DoubleReal > intensity
Peak intensity (0 if peak index is -1 or -2)
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:307
Size getPeakCount() const
Returns the peak count of all traces.
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:163
void resize(Size s)
Definition: MSExperiment.h:122
virtual DoubleReal computeTheoretical(const FeatureFinderAlgorithmPickedHelperStructs::MassTrace< PeakType > &trace, Size k)=0
void setIsotopeDistribution(const IsotopeDistribution &isotopes)
sets the isotope distribution of the element
CoordinateType getMaxMZ() const
returns the maximal m/z value
Definition: MSExperiment.h:527
std::vector< DoubleReal > theoretical_mz
Theoretical m/z value of the isotope peak.
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:311
void setParameters(const Param ¶m)
Sets the parameters.
DoubleReal intensity_rt_step_
RT bin width.
Definition: FeatureFinderAlgorithmPicked.h:1150
ConstAreaIterator areaEndConst() const
Returns an non-mutable invalid area iterator marking the end of an area.
Definition: MSExperiment.h:337
#define LOG_DEBUG
Macro for general debugging information.
Definition: LogStream.h:459
The purpose of this struct is to provide definitions of classes and typedefs which are used throughou...
Definition: FeatureFinderDefs.h:51
CoordinateType getMinMZ() const
returns the minimal m/z value
Definition: MSExperiment.h:521
SpectrumType::FloatDataArrays FloatDataArrays
Definition: FeatureFinderAlgorithmPicked.h:94
CoordinateType getMinRT() const
returns the minimal retention time value
Definition: MSExperiment.h:533
const double PROTON_MASS_U
void setIntensity(IntensityType intensity)
Non-mutable access to the data point intensity (height)
Definition: Peak2D.h:167
std::vector< FloatDataArray > FloatDataArrays
Float data array vector type.
Definition: MSSpectrum.h:113
TheoreticalIsotopePattern theoretical_pattern
Theoretical isotope pattern.
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:313
void store(const String &filename)
Writes the data to a file.
static const String getProductName()
Definition: FeatureFinderAlgorithmPicked.h:1108
Fitter for RT profiles using a gaussian background model.
Definition: GaussTraceFitter.h:53
std::vector< Size > spectrum
Spectrum index (undefined if peak index is -1 or -2)
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:305
Size nearest_(DoubleReal pos, const SpectrumType &spec, Size start) const
Returns the index of the peak nearest to m/z pos in spectrum spec (linear search starting from index ...
Definition: FeatureFinderAlgorithmPicked.h:1612
bool isValid() const
Checks if this Trace is valid (has more than 2 points)
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:142
File adapter for MzML files.
Definition: MzMLFile.h:58
Representation of an element.
Definition: Element.h:54
UInt intensity_bins_
Number of bins (in RT and MZ) for intensity significance estimation.
Definition: FeatureFinderAlgorithmPicked.h:1138
DoubleReal max_rt_span_
Maximum RT range the model is allowed to span.
Definition: FeatureFinderAlgorithmPicked.h:1142
void writeFeatureDebugInfo_(TraceFitter< PeakType > *fitter, const MassTraces &traces, const MassTraces &new_traces, bool feature_ok, const String error_msg, const DoubleReal final_score, const Int plot_nr, const PeakType &peak, const String path="debug/features/")
Creates several files containing plots and viewable data of the fitted mass trace.
Definition: FeatureFinderAlgorithmPicked.h:2181
void endProgress() const
Ends the progress display.
FeatureFinderAlgorithmPickedHelperStructs::IsotopePattern IsotopePattern
Definition: FeatureFinderAlgorithmPicked.h:107
CoordinateType width() const
Returns the width of the area i.e. the difference of dimension zero (X).
Definition: DIntervalBase.h:292
const DataValue & getValue(const String &key) const
Returns a value of a parameter.
void setRT(CoordinateType coordinate)
Mutable access to the RT coordinate (index 0)
Definition: Peak2D.h:209
static DoubleReal pearsonCorrelationCoefficient(IteratorType1 begin_a, IteratorType1 end_a, IteratorType2 begin_b, IteratorType2 end_b)
Calculates the Pearson correlation coefficient for the values in [begin_a, end_a) and [begin_b...
Definition: StatisticFunctions.h:285
ConstAreaIterator areaBeginConst(CoordinateType min_rt, CoordinateType max_rt, CoordinateType min_mz, CoordinateType max_mz) const
Returns a non-mutable area iterator for area.
Definition: MSExperiment.h:327
DoubleReal theoretical_int
Theoretical intensity value (scaled to [0,1])
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:91
const std::vector< Feature > & getSubordinates() const
immutable access to subordinate features
void setMaxIsotope(Size max_isotope)
sets the maximal isotope with max_isotope
const Element * getElement(const String &name) const
FeatureFinderAlgorithmPicked()
default constructor
Definition: FeatureFinderAlgorithmPicked.h:111
DoubleReal optional_fit_improvement_
Minimal imrovment for leaving out optional isotope.
Definition: FeatureFinderAlgorithmPicked.h:1136
const std::vector< ConvexHull2D > & getConvexHulls() const
Non-mutable access to the convex hulls.
virtual bool checkMaximalRTSpan(const DoubleReal max_rt_span)=0
bool checkFeatureQuality_(TraceFitter< PeakType > *fitter, MassTraces &feature_traces, const DoubleReal &seed_mz, const DoubleReal &min_feature_score, String &error_msg, DoubleReal &fit_score, DoubleReal &correlation, DoubleReal &final_score)
Checks the feature based on different score thresholds and model constraints.
Definition: FeatureFinderAlgorithmPicked.h:2080
Iterator begin()
Definition: IsotopeDistribution.h:185
virtual DoubleReal getHeight() const =0
virtual DoubleReal getFWHM() const =0
void store(const String &filename, const MapType &map) const
Stores a map in a MzML file.
Definition: MzMLFile.h:126
void trimLeft(DoubleReal cutoff)
Trims the left side of the isotope distribution to isotopes with a significant contribution.
DoubleReal max_feature_intersection_
Maximum allowed feature intersection (if larger, that one of the feature is removed) ...
Definition: FeatureFinderAlgorithmPicked.h:1143
DoubleReal min_trace_score_
Minimum quality of a traces.
Definition: FeatureFinderAlgorithmPicked.h:1140
virtual DoubleReal getUpperRTBound() const =0
static StringList create(const String &list, const char splitter= ',')
Returns a list that is created by splitting the given (comma-separated) string (String are not trimme...
void estimateFromPeptideWeight(double average_weight)
Estimate Peptide Isotopedistribution from weight and number of isotopes that should be reported...
void trimRight(DoubleReal cutoff)
Trims the right side of the isotope distribution to isotopes with a significant contribution.
DBoundingBox< 2 > getBoundingBox() const
returns the bounding box of the feature hull points
DoubleReal findBestIsotopeFit_(const Seed ¢er, UInt charge, IsotopePattern &best_pattern) const
Finds the best fitting position of the isotopic pattern estimate defined by center.
Definition: FeatureFinderAlgorithmPicked.h:1270
DoubleReal baseline
Estimated baseline in the region of the feature (used for the fit)
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:271
An LC-MS feature.
Definition: Feature.h:66
FeatureFinderAlgorithm< PeakType, FeatureType >::MapType MapType
Definition: FeatureFinderAlgorithmPicked.h:91
Helper structure for a found isotope pattern used in FeatureFinderAlgorithmPicked.
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:300
FeatureFinderAlgorithmPickedHelperStructs::Seed Seed
Definition: FeatureFinderAlgorithmPicked.h:103
void setWidth(WidthType fwhm)
Set the width of the feature (FWHM)
DoubleReal isotopeScore_(const TheoreticalIsotopePattern &isotopes, IsotopePattern &pattern, bool consider_mz_distances) const
Calculates a score between 0 and 1 for the correlation between theoretical and found isotope pattern...
Definition: FeatureFinderAlgorithmPicked.h:1738
void setValidStrings(const String &key, const std::vector< String > &strings)
Sets the valid strings for the parameter key.
virtual void fit(FeatureFinderAlgorithmPickedHelperStructs::MassTraces< PeakType > &traces)=0
virtual void run()
Main method for actual FeatureFinder.
Definition: FeatureFinderAlgorithmPicked.h:218
Management and storage of parameters / INI files.
Definition: Param.h:69
CoordinateType getMZ() const
Returns the m/z coordinate (index 1)
Definition: Peak2D.h:191
Invalid value exception.
Definition: Exception.h:336
DoubleReal intensity_percentage_optional_
Isotope pattern intensity contribution of optional peaks.
Definition: FeatureFinderAlgorithmPicked.h:1135
Representation of a mass spectrometry experiment.
Definition: MSExperiment.h:68
void swapFeaturesOnly(FeatureMap &from)
Swaps the feature content (plus its range information) of this map with the content of from...
Definition: FeatureMap.h:377
UInt max_missing_trace_peaks_
Stores mass_trace:max_missing.
Definition: FeatureFinderAlgorithmPicked.h:1132
std::map< String, UInt > aborts_
Array of abort reasons.
Definition: FeatureFinderAlgorithmPicked.h:1121
DoubleReal slope_bound_
Max slope of mass trace intensities.
Definition: FeatureFinderAlgorithmPicked.h:1133
PositionType const & minPosition() const
Accessor to minimum position.
Definition: DIntervalBase.h:121
FeatureFinderAlgorithmPickedHelperStructs::MassTraces< PeakType > MassTraces
Definition: FeatureFinderAlgorithmPicked.h:105
This class provides Input/Output functionality for feature maps.
Definition: FeatureXMLFile.h:59
void setOverallQuality(QualityType q)
Set the overall quality.
std::vector< TheoreticalIsotopePattern > isotope_distributions_
Vector of precalculated isotope distributions for several mass windows.
Definition: FeatureFinderAlgorithmPicked.h:1158
static const ElementDB * getInstance()
returns a pointer to the singleton instance of the element db
Definition: ElementDB.h:78
CoordinateType getMaxRT() const
returns the maximal retention time value
Definition: MSExperiment.h:539
void setMinInt(const String &key, Int min)
Sets the minimum value for the integer or integer list parameter key.
size_t Size
Size type e.g. used as variable which can hold result of size()
Definition: Types.h:144
std::vector< std::pair< DoubleReal, const PeakType * > > peaks
Contained peaks (pair of RT and pointer to peak)
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:94
UInt min_spectra_
Number of spectra that have to show the same mass (for finding a mass trace)
Definition: FeatureFinderAlgorithmPicked.h:1131
Size optional_end
Number of optional peaks at the end of the pattern.
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:284
DoubleReal intensityScore_(Size rt_bin, Size mz_bin, DoubleReal intensity) const
Definition: FeatureFinderAlgorithmPicked.h:1930
void startProgress(SignedSize begin, SignedSize end, const String &label) const
Initializes the progress display.
bool empty() const
Definition: MSExperiment.h:127
void extendMassTrace_(MassTrace &trace, SignedSize spectrum_index, DoubleReal mz, bool increase_rt, Size meta_index_overall, DoubleReal min_rt=0.0, DoubleReal max_rt=0.0) const
Extends a single mass trace in one RT direction.
Definition: FeatureFinderAlgorithmPicked.h:1510
MapType map_
editable copy of the map
Definition: FeatureFinderAlgorithmPicked.h:1115
void setCharge(const ChargeType &ch)
Set charge state.
Size trimmed_left
The number of isotopes trimmed on the left side. This is needed to reconstruct the monoisotopic peak...
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:288
FeatureFinderAlgorithmPickedHelperStructs::MassTrace< PeakType > MassTrace
Definition: FeatureFinderAlgorithmPicked.h:104
void setProgress(SignedSize value) const
Sets the current progress.
PositionType const & maxPosition() const
Accessor to maximum position.
Definition: DIntervalBase.h:127
Size getTheoreticalmaxPosition() const
Returns the theoretical maximum trace index.
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:195
A RT Profile fitter using an Exponential Gaussian Hybrid background model.
Definition: EGHTraceFitter.h:58
virtual void updateMembers_()
This method is used to update extra member variables at the end of the setParameters() method...
Definition: FeatureFinderAlgorithmPicked.h:1161
A D-dimensional bounding box.
Definition: DBoundingBox.h:51
DoubleReal getRT() const
Definition: MSSpectrum.h:215
FeatureMapType seeds_
User-specified seed list.
Definition: FeatureFinderAlgorithmPicked.h:1125
void updateBaseline()
Sets the baseline to the lowest contained peak of the trace.
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:216
virtual DoubleReal getCenter() const =0
DoubleReal intensity_mz_step_
m/z bin width
Definition: FeatureFinderAlgorithmPicked.h:1152
int Int
Signed integer type.
Definition: Types.h:100
void set(const ContainerType &distribution)
overwrites the container which holds the distribution using distribution
void setMinFloat(const String &key, DoubleReal min)
Sets the minimum value for the floating point or floating point list parameter key.
static FeatureFinderAlgorithm< PeakType, FeatureType > * create()
Definition: FeatureFinderAlgorithmPicked.h:1103
virtual void setSeeds(const FeatureMapType &seeds)
Sets a reference to the calling FeatureFinder.
Definition: FeatureFinderAlgorithmPicked.h:212
std::vector< DoubleReal > intensity
Vector of intensity contributions.
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:280
void abort_(const Seed &seed, const String &reason)
Writes the abort reason to the log file and counts occurrences for each reason.
Definition: FeatureFinderAlgorithmPicked.h:1182
void sortByIntensity(bool reverse=false)
Sorts the peaks according to ascending intensity.
Definition: FeatureMap.h:297
std::pair< DoubleReal, DoubleReal > getRTBounds() const
Returns the RT boundaries of the mass traces.
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:246
Size peak
Peak index.
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:67
std::vector< SignedSize > peak
Peak index (-1 if peak was not found, -2 if it was removed to improve the isotope fit) ...
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:303
This class provides some basic file handling methods for text files.
Definition: TextFile.h:47
FeatureFinderAlgorithmPickedHelperStructs::TheoreticalIsotopePattern TheoreticalIsotopePattern
Definition: FeatureFinderAlgorithmPicked.h:106
void defaultsToParam_()
Updates the parameters after the defaults have been set in the constructor.
Real intensity
Intensity.
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:69
FeatureFinderAlgorithmPicked & operator=(const FeatureFinderAlgorithmPicked &)
Not implemented.
void setMaxFloat(const String &key, DoubleReal max)
Sets the maximum value for the floating point or floating point list parameter key.
Helper struct for mass traces used in FeatureFinderAlgorithmPicked.
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:83
Helper structure for seeds used in FeatureFinderAlgorithmPicked.
Definition: FeatureFinderAlgorithmPickedHelperStructs.h:62
double DoubleReal
Double-precision real type.
Definition: Types.h:118