00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035 #ifndef OPENMS_TRANSFORMATIONS_FEATUREFINDER_FEATUREFINDERALGORITHMPICKED_H
00036 #define OPENMS_TRANSFORMATIONS_FEATUREFINDER_FEATUREFINDERALGORITHMPICKED_H
00037
00038 #include <OpenMS/TRANSFORMATIONS/FEATUREFINDER/FeatureFinderAlgorithm.h>
00039 #include <OpenMS/TRANSFORMATIONS/FEATUREFINDER/FeatureFinderAlgorithmPickedHelperStructs.h>
00040 #include <OpenMS/TRANSFORMATIONS/FEATUREFINDER/TraceFitter.h>
00041 #include <OpenMS/TRANSFORMATIONS/FEATUREFINDER/EGHTraceFitter.h>
00042 #include <OpenMS/TRANSFORMATIONS/FEATUREFINDER/GaussTraceFitter.h>
00043
00044 #include <OpenMS/FORMAT/MzMLFile.h>
00045 #include <OpenMS/FORMAT/FeatureXMLFile.h>
00046 #include <OpenMS/FORMAT/TextFile.h>
00047 #include <OpenMS/CHEMISTRY/IsotopeDistribution.h>
00048 #include <OpenMS/MATH/STATISTICS/StatisticFunctions.h>
00049 #include <OpenMS/MATH/MISC/MathFunctions.h>
00050 #include <OpenMS/CONCEPT/Constants.h>
00051 #include <OpenMS/CHEMISTRY/Element.h>
00052 #include <OpenMS/CHEMISTRY/ElementDB.h>
00053 #include <OpenMS/CHEMISTRY/IsotopeDistribution.h>
00054
00055 #include <boost/math/special_functions/fpclassify.hpp>
00056
00057 #include <numeric>
00058 #include <fstream>
00059 #include <algorithm>
00060
00061 #include <gsl/gsl_rng.h>
00062 #include <gsl/gsl_vector.h>
00063 #include <gsl/gsl_multifit_nlin.h>
00064 #include <gsl/gsl_blas.h>
00065
00066 #include <QtCore/QDir>
00067
00068 #ifdef _OPENMP
00069 #include <omp.h>
00070 #endif
00071
00072
00073 namespace OpenMS
00074 {
00088 template <class PeakType, class FeatureType>
00089 class FeatureFinderAlgorithmPicked :
00090 public FeatureFinderAlgorithm<PeakType, FeatureType>,
00091 public FeatureFinderDefs
00092 {
00093 public:
00095
00096 typedef typename FeatureFinderAlgorithm<PeakType, FeatureType>::MapType MapType;
00097 typedef typename FeatureFinderAlgorithm<PeakType, FeatureType>::FeatureMapType FeatureMapType;
00098 typedef typename MapType::SpectrumType SpectrumType;
00099 typedef typename SpectrumType::FloatDataArrays FloatDataArrays;
00101
00102 using FeatureFinderAlgorithm<PeakType, FeatureType>::param_;
00103 using FeatureFinderAlgorithm<PeakType, FeatureType>::features_;
00104 using FeatureFinderAlgorithm<PeakType, FeatureType>::ff_;
00105 using FeatureFinderAlgorithm<PeakType, FeatureType>::defaults_;
00106
00107 protected:
00108 typedef FeatureFinderAlgorithmPickedHelperStructs::Seed Seed;
00109 typedef typename FeatureFinderAlgorithmPickedHelperStructs::MassTrace<PeakType> MassTrace;
00110 typedef typename FeatureFinderAlgorithmPickedHelperStructs::MassTraces<PeakType> MassTraces;
00111 typedef FeatureFinderAlgorithmPickedHelperStructs::TheoreticalIsotopePattern TheoreticalIsotopePattern;
00112 typedef FeatureFinderAlgorithmPickedHelperStructs::IsotopePattern IsotopePattern;
00113
00114 public:
00116 FeatureFinderAlgorithmPicked() :
00117 FeatureFinderAlgorithm<PeakType, FeatureType>(),
00118 map_(),
00119 log_()
00120 {
00121
00122 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).");
00123 defaults_.setValidStrings("debug", StringList::create("true,false"));
00124
00125 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.");
00126 defaults_.setMinInt("intensity:bins", 1);
00127 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)");
00128
00129 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!");
00130 defaults_.setMinFloat("mass_trace:mz_tolerance", 0.0);
00131 defaults_.setValue("mass_trace:min_spectra", 10, "Number of spectra that have to show a similar peak mass in a mass trace.");
00132 defaults_.setMinInt("mass_trace:min_spectra", 1);
00133 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'!");
00134 defaults_.setMinInt("mass_trace:max_missing", 0);
00135 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.");
00136 defaults_.setMinFloat("mass_trace:slope_bound", 0.0);
00137 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).");
00138
00139 defaults_.setValue("isotopic_pattern:charge_low", 1, "Lowest charge to search for.");
00140 defaults_.setMinInt("isotopic_pattern:charge_low", 1);
00141 defaults_.setValue("isotopic_pattern:charge_high", 4, "Highest charge to search for.");
00142 defaults_.setMinInt("isotopic_pattern:charge_high", 1);
00143 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!");
00144 defaults_.setMinFloat("isotopic_pattern:mz_tolerance", 0.0);
00145 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"));
00146 defaults_.setMinFloat("isotopic_pattern:intensity_percentage", 0.0);
00147 defaults_.setMaxFloat("isotopic_pattern:intensity_percentage", 100.0);
00148 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"));
00149 defaults_.setMinFloat("isotopic_pattern:intensity_percentage_optional", 0.0);
00150 defaults_.setMaxFloat("isotopic_pattern:intensity_percentage_optional", 100.0);
00151 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"));
00152 defaults_.setMinFloat("isotopic_pattern:optional_fit_improvement", 0.0);
00153 defaults_.setMaxFloat("isotopic_pattern:optional_fit_improvement", 100.0);
00154 defaults_.setValue("isotopic_pattern:mass_window_width", 25.0, "Window width in Dalton for precalculation of estimated isotope distributions.", StringList::create("advanced"));
00155 defaults_.setMinFloat("isotopic_pattern:mass_window_width", 1.0);
00156 defaults_.setMaxFloat("isotopic_pattern:mass_window_width", 200.0);
00157 defaults_.setValue("isotopic_pattern:abundance_12C", 98.93, "Rel. abundance of the light carbon. Modify if labeled.", StringList::create("advanced"));
00158 defaults_.setMinFloat("isotopic_pattern:abundance_12C", 0.0);
00159 defaults_.setMaxFloat("isotopic_pattern:abundance_12C", 100.0);
00160 defaults_.setValue("isotopic_pattern:abundance_14N", 99.632, "Rel. abundance of the light nitrogen. Modify if labeled.", StringList::create("advanced"));
00161 defaults_.setMinFloat("isotopic_pattern:abundance_14N", 0.0);
00162 defaults_.setMaxFloat("isotopic_pattern:abundance_14N", 100.0);
00163
00164 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).");
00165
00166 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.");
00167 defaults_.setMinFloat("seed:min_score", 0.0);
00168 defaults_.setMaxFloat("seed:min_score", 1.0);
00169 defaults_.setSectionDescription("seed", "Settings that determine which peaks are considered a seed");
00170
00171 defaults_.setValue("fit:epsilon_abs", 0.0001, "Absolute epsilon used for convergence of the fit.", StringList::create("advanced"));
00172 defaults_.setMinFloat("fit:epsilon_abs", 0.0);
00173 defaults_.setValue("fit:epsilon_rel", 0.0001, "Relative epsilon used for convergence of the fit.", StringList::create("advanced"));
00174 defaults_.setMinFloat("fit:epsilon_rel", 0.0);
00175 defaults_.setValue("fit:max_iterations", 500, "Maximum number of iterations of the fit.", StringList::create("advanced"));
00176 defaults_.setMinInt("fit:max_iterations", 1);
00177 defaults_.setSectionDescription("fit", "Settings for the model fitting");
00178
00179 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.");
00180 defaults_.setMinFloat("feature:min_score", 0.0);
00181 defaults_.setMaxFloat("feature:min_score", 1.0);
00182 defaults_.setValue("feature:min_isotope_fit", 0.8, "Minimum isotope fit of the feature before model fitting.", StringList::create("advanced"));
00183 defaults_.setMinFloat("feature:min_isotope_fit", 0.0);
00184 defaults_.setMaxFloat("feature:min_isotope_fit", 1.0);
00185 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"));
00186 defaults_.setMinFloat("feature:min_trace_score", 0.0);
00187 defaults_.setMaxFloat("feature:min_trace_score", 1.0);
00188 defaults_.setValue("feature:min_rt_span", 0.333, "Minimum RT span in relation to extended area that has to remain after model fitting.", StringList::create("advanced"));
00189 defaults_.setMinFloat("feature:min_rt_span", 0.0);
00190 defaults_.setMaxFloat("feature:min_rt_span", 1.0);
00191 defaults_.setValue("feature:max_rt_span", 2.5, "Maximum RT span in relation to extended area that the model is allowed to have.", StringList::create("advanced"));
00192 defaults_.setMinFloat("feature:max_rt_span", 0.5);
00193 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"));
00194 defaults_.setValidStrings("feature:rt_shape", StringList::create("symmetric,asymmetric"));
00195 defaults_.setValue("feature:max_intersection", 0.35, "Maximum allowed intersection of features.", StringList::create("advanced"));
00196 defaults_.setMinFloat("feature:max_intersection", 0.0);
00197 defaults_.setMaxFloat("feature:max_intersection", 1.0);
00198 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.");
00199 defaults_.setValidStrings("feature:reported_mz", StringList::create("maximum,average,monoisotopic"));
00200 defaults_.setSectionDescription("feature", "Settings for the features (intensity, quality assessment, ...)");
00201
00202 defaults_.setValue("user-seed:rt_tolerance", 5.0, "Allowed RT deviation of seeds from the user-specified seed position.");
00203 defaults_.setMinFloat("user-seed:rt_tolerance", 0.0);
00204 defaults_.setValue("user-seed:mz_tolerance", 1.1, "Allowed m/z deviation of seeds from the user-specified seed position.");
00205 defaults_.setMinFloat("user-seed:mz_tolerance", 0.0);
00206 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.");
00207 defaults_.setMinFloat("user-seed:min_score", 0.0);
00208 defaults_.setMaxFloat("user-seed:min_score", 1.0);
00209 defaults_.setSectionDescription("user-seed", "Settings for user-specified seeds.");
00210
00211 defaults_.setValue("debug:pseudo_rt_shift", 500.0, "Pseudo RT shift used when .", StringList::create("advanced"));
00212 defaults_.setMinFloat("debug:pseudo_rt_shift", 1.0);
00213 this->defaultsToParam_();
00214 }
00215
00216
00217 virtual void setSeeds(const FeatureMapType& seeds)
00218 {
00219 seeds_ = seeds;
00220 }
00221
00223 virtual void run()
00224 {
00225
00226
00227
00228
00229
00230 DoubleReal min_feature_score = param_.getValue("feature:min_score");
00231
00232 SignedSize charge_low = (Int)param_.getValue("isotopic_pattern:charge_low");
00233 SignedSize charge_high = (Int)param_.getValue("isotopic_pattern:charge_high");
00234
00235 UInt max_iterations = param_.getValue("fit:max_iterations");
00236 DoubleReal epsilon_abs = param_.getValue("fit:epsilon_abs");
00237 DoubleReal epsilon_rel = param_.getValue("fit:epsilon_rel");
00238
00239 Size max_isotopes = 20;
00240
00241
00242 DoubleReal abundance_12C = param_.getValue("isotopic_pattern:abundance_12C");
00243 DoubleReal abundance_14N = param_.getValue("isotopic_pattern:abundance_14N");
00244
00245 const Element * carbon_const = ElementDB::getInstance()->getElement("Carbon");
00246 Element * carbon = const_cast<Element *>(carbon_const);
00247
00248 if (param_.getValue("isotopic_pattern:abundance_12C") != defaults_.getValue("isotopic_pattern:abundance_12C"))
00249 {
00250 max_isotopes += 1000;
00251 IsotopeDistribution isotopes;
00252 std::vector<std::pair<Size, double> > container;
00253 container.push_back(std::make_pair(12, abundance_12C / 100.0));
00254 container.push_back(std::make_pair(13, 1.0 - (abundance_12C / 100.0)));
00255 isotopes.set(container);
00256 carbon->setIsotopeDistribution(isotopes);
00257 }
00258
00259 const Element * nitrogen_const = ElementDB::getInstance()->getElement("Nitrogen");
00260 Element * nitrogen = const_cast<Element *>(nitrogen_const);
00261
00262 if (param_.getValue("isotopic_pattern:abundance_14N") != defaults_.getValue("isotopic_pattern:abundance_14N"))
00263 {
00264 max_isotopes += 1000;
00265 IsotopeDistribution isotopes;
00266 std::vector<std::pair<Size, double> > container;
00267 container.push_back(std::make_pair(14, abundance_14N / 100.0));
00268 container.push_back(std::make_pair(15, 1.0 - (abundance_14N / 100.0)));
00269 isotopes.set(container);
00270 nitrogen->setIsotopeDistribution(isotopes);
00271 }
00272
00273
00274
00275 Param trace_fitter_params;
00276 trace_fitter_params.setValue("max_iteration", max_iterations);
00277 trace_fitter_params.setValue("epsilon_abs", epsilon_abs);
00278 trace_fitter_params.setValue("epsilon_rel", epsilon_rel);
00279
00280
00281 map_ = *(FeatureFinderAlgorithm<PeakType, FeatureType>::map_);
00282
00283
00284 bool user_seeds = (seeds_.size() > 0);
00285 if (user_seeds)
00286 {
00287 seeds_.sortByMZ();
00288 }
00289 DoubleReal user_rt_tol = param_.getValue("user-seed:rt_tolerance");
00290 DoubleReal user_mz_tol = param_.getValue("user-seed:mz_tolerance");
00291 DoubleReal user_seed_score = param_.getValue("user-seed:min_score");
00292
00293
00294 UInt charge_count = charge_high - charge_low + 1;
00295 for (Size s = 0; s < map_.size(); ++s)
00296 {
00297 Size scan_size = map_[s].size();
00298 map_[s].getFloatDataArrays().resize(3 + 2 * charge_count);
00299 map_[s].getFloatDataArrays()[0].setName("trace_score");
00300 map_[s].getFloatDataArrays()[0].assign(scan_size, 0.0);
00301 map_[s].getFloatDataArrays()[1].setName("intensity_score");
00302 map_[s].getFloatDataArrays()[1].assign(scan_size, 0.0);
00303 map_[s].getFloatDataArrays()[2].setName("local_max");
00304 map_[s].getFloatDataArrays()[2].assign(scan_size, 0.0);
00305
00306 UInt charge = charge_low;
00307 for (Size i = 3; i < 3 + charge_count; ++i)
00308 {
00309 map_[s].getFloatDataArrays()[i].setName(String("pattern_score_") + charge);
00310 map_[s].getFloatDataArrays()[i].assign(scan_size, 0.0);
00311 ++charge;
00312 }
00313
00314 charge = charge_low;
00315 for (Size i = 3 + charge_count; i < 3 + 2 * charge_count; ++i)
00316 {
00317 map_[s].getFloatDataArrays()[i].setName(String("overall_score_") + charge);
00318 map_[s].getFloatDataArrays()[i].assign(scan_size, 0.0);
00319 ++charge;
00320 }
00321 }
00322
00323 int gl_progress = 0;
00324 debug_ = ((String)(param_.getValue("debug")) == "true");
00325
00326 if (debug_)
00327 {
00328 QDir dir(".");
00329 dir.mkpath("debug/features");
00330 log_.open("debug/log.txt");
00331 }
00332
00333
00334
00335
00336
00337 if (debug_) log_ << "Precalculating intensity thresholds ..." << std::endl;
00338
00339 {
00340 ff_->startProgress(0, intensity_bins_ * intensity_bins_, "Precalculating intensity scores");
00341 DoubleReal rt_start = map_.getMinRT();
00342 DoubleReal mz_start = map_.getMinMZ();
00343 intensity_rt_step_ = (map_.getMaxRT() - rt_start) / (DoubleReal)intensity_bins_;
00344 intensity_mz_step_ = (map_.getMaxMZ() - mz_start) / (DoubleReal)intensity_bins_;
00345 intensity_thresholds_.resize(intensity_bins_);
00346 for (Size rt = 0; rt < intensity_bins_; ++rt)
00347 {
00348 intensity_thresholds_[rt].resize(intensity_bins_);
00349 DoubleReal min_rt = rt_start + rt * intensity_rt_step_;
00350 DoubleReal max_rt = rt_start + (rt + 1) * intensity_rt_step_;
00351 std::vector<DoubleReal> tmp;
00352 for (Size mz = 0; mz < intensity_bins_; ++mz)
00353 {
00354 ff_->setProgress(rt * intensity_bins_ + mz);
00355 DoubleReal min_mz = mz_start + mz * intensity_mz_step_;
00356 DoubleReal max_mz = mz_start + (mz + 1) * intensity_mz_step_;
00357
00358
00359 tmp.clear();
00360 for (typename MapType::ConstAreaIterator it = map_.areaBeginConst(min_rt, max_rt, min_mz, max_mz); it != map_.areaEndConst(); ++it)
00361 {
00362 tmp.push_back(it->getIntensity());
00363 }
00364
00365 intensity_thresholds_[rt][mz].assign(21, 0.0);
00366
00367 if (!tmp.empty())
00368 {
00369 std::sort(tmp.begin(), tmp.end());
00370 for (Size i = 0; i < 21; ++i)
00371 {
00372 Size index = (Size) std::floor(0.05 * i * (tmp.size() - 1));
00373 intensity_thresholds_[rt][mz][i] = tmp[index];
00374 }
00375 }
00376 }
00377 }
00378
00379
00380 for (Size s = 0; s < map_.size(); ++s)
00381 {
00382 for (Size p = 0; p < map_[s].size(); ++p)
00383 {
00384 map_[s].getFloatDataArrays()[1][p] = intensityScore_(s, p);
00385 }
00386 }
00387 ff_->endProgress();
00388 }
00389
00390
00391
00392
00393
00394
00395 {
00396 Size end_iteration = map_.size() - std::min((Size) min_spectra_, map_.size());
00397 ff_->startProgress(min_spectra_, end_iteration, "Precalculating mass trace scores");
00398
00399 for (Size s = min_spectra_; s < end_iteration; ++s)
00400 {
00401 ff_->setProgress(s);
00402 const SpectrumType& spectrum = map_[s];
00403
00404 for (Size p = 0; p < spectrum.size(); ++p)
00405 {
00406 std::vector<DoubleReal> scores;
00407 scores.reserve(2 * min_spectra_);
00408
00409 DoubleReal pos = spectrum[p].getMZ();
00410 Real inte = spectrum[p].getIntensity();
00411
00412
00413 bool is_max_peak = true;
00414 for (Size i = 1; i <= min_spectra_; ++i)
00415 {
00416 try
00417 {
00418 Size spec_index = map_[s + i].findNearest(pos);
00419 DoubleReal position_score = positionScore_(pos, map_[s + i][spec_index].getMZ(), trace_tolerance_);
00420 if (position_score > 0 && map_[s + i][spec_index].getIntensity() > inte) is_max_peak = false;
00421 scores.push_back(position_score);
00422 }
00423 catch (...)
00424 {
00425 scores.push_back(0.0);
00426 }
00427 }
00428 for (Size i = 1; i <= min_spectra_; ++i)
00429 {
00430 try
00431 {
00432 Size spec_index = map_[s - i].findNearest(pos);
00433 DoubleReal position_score = positionScore_(pos, map_[s - i][spec_index].getMZ(), trace_tolerance_);
00434 if (position_score > 0 && map_[s - i][spec_index].getIntensity() > inte) is_max_peak = false;
00435 scores.push_back(position_score);
00436 }
00437 catch (...)
00438 {
00439 scores.push_back(0.0);
00440 }
00441 }
00442
00443 DoubleReal trace_score = std::accumulate(scores.begin(), scores.end(), 0.0) / scores.size();
00444
00445
00446 map_[s].getFloatDataArrays()[0][p] = trace_score;
00447 map_[s].getFloatDataArrays()[2][p] = is_max_peak;
00448 }
00449 }
00450 ff_->endProgress();
00451 }
00452
00453
00454
00455
00456
00457
00458 {
00459 DoubleReal max_mass = map_.getMaxMZ() * charge_high;
00460 Size num_isotopes = std::ceil(max_mass / mass_window_width_) + 1;
00461 ff_->startProgress(0, num_isotopes, "Precalculating isotope distributions");
00462
00463
00464 isotope_distributions_.resize(num_isotopes);
00465
00466
00467 for (Size index = 0; index < num_isotopes; ++index)
00468 {
00469
00470 IsotopeDistribution d;
00471 d.setMaxIsotope(max_isotopes);
00472 d.estimateFromPeptideWeight(0.5 * mass_window_width_ + index * mass_window_width_);
00473
00474 Size size_before = d.size();
00475 d.trimLeft(intensity_percentage_optional_);
00476 isotope_distributions_[index].trimmed_left = size_before - d.size();
00477 d.trimRight(intensity_percentage_optional_);
00478
00479 for (IsotopeDistribution::Iterator it = d.begin(); it != d.end(); ++it)
00480 {
00481 isotope_distributions_[index].intensity.push_back(it->second);
00482
00483 }
00484
00485
00486 Size begin = 0;
00487 Size end = 0;
00488 bool is_begin = true;
00489 bool is_end = false;
00490
00491 for (Size i = 0; i < isotope_distributions_[index].intensity.size(); ++i)
00492 {
00493 if (isotope_distributions_[index].intensity[i] < intensity_percentage_)
00494 {
00495 if (!is_end && !is_begin) is_end = true;
00496 if (is_begin) ++begin;
00497 else if (is_end) ++end;
00498 }
00499 else if (is_begin)
00500 {
00501 is_begin = false;
00502 }
00503 }
00504 isotope_distributions_[index].optional_begin = begin;
00505 isotope_distributions_[index].optional_end = end;
00506
00507 DoubleReal max = 0.0;
00508 for (Size i = 0; i < isotope_distributions_[index].intensity.size(); ++i)
00509 {
00510 if (isotope_distributions_[index].intensity[i] > max)
00511 {
00512 max = isotope_distributions_[index].intensity[i];
00513 }
00514 }
00515 isotope_distributions_[index].max = max;
00516 for (Size i = 0; i < isotope_distributions_[index].intensity.size(); ++i)
00517 {
00518 isotope_distributions_[index].intensity[i] /= max;
00519 }
00520
00521
00522 }
00523
00524 ff_->endProgress();
00525 }
00526
00527
00528
00529
00530
00531 Int plot_nr_global = -1;
00532 Int feature_nr_global = 0;
00533 for (SignedSize c = charge_low; c <= charge_high; ++c)
00534 {
00535 UInt meta_index_isotope = 3 + c - charge_low;
00536 UInt meta_index_overall = 3 + charge_count + c - charge_low;
00537
00538 Size feature_candidates = 0;
00539 std::vector<Seed> seeds;
00540
00541
00542
00543
00544 ff_->startProgress(0, map_.size(), String("Calculating isotope pattern scores for charge ") + String(c));
00545 for (Size s = 0; s < map_.size(); ++s)
00546 {
00547 ff_->setProgress(s);
00548 const SpectrumType& spectrum = map_[s];
00549 for (Size p = 0; p < spectrum.size(); ++p)
00550 {
00551 DoubleReal mz = spectrum[p].getMZ();
00552
00553
00554 const TheoreticalIsotopePattern& isotopes = getIsotopeDistribution_(mz * c);
00555
00556 Size max_isotope = std::max_element(isotopes.intensity.begin(), isotopes.intensity.end()) - isotopes.intensity.begin();
00557
00558 Size peak_index = spectrum.findNearest(mz - ((DoubleReal)(isotopes.size() + 1) / c));
00559 IsotopePattern pattern(isotopes.size());
00560
00561 for (Size i = 0; i < isotopes.size(); ++i)
00562 {
00563 DoubleReal isotope_pos = mz + ((DoubleReal)i - max_isotope) / c;
00564 findIsotope_(isotope_pos, s, pattern, i, peak_index);
00565 }
00566
00567 DoubleReal pattern_score = isotopeScore_(isotopes, pattern, true);
00568
00569
00570 if (pattern_score > 0.0)
00571 {
00572 for (Size i = 0; i < pattern.peak.size(); ++i)
00573 {
00574 if (pattern.peak[i] >= 0 && pattern_score > map_[pattern.spectrum[i]].getFloatDataArrays()[meta_index_isotope][pattern.peak[i]])
00575 {
00576 map_[pattern.spectrum[i]].getFloatDataArrays()[meta_index_isotope][pattern.peak[i]] = pattern_score;
00577 }
00578 }
00579 }
00580 }
00581 }
00582 ff_->endProgress();
00583
00584
00585
00586
00587 Size end_of_iteration = map_.size() - std::min((Size) min_spectra_, map_.size());
00588 ff_->startProgress(min_spectra_, end_of_iteration, String("Finding seeds for charge ") + String(c));
00589
00590 DoubleReal min_seed_score = param_.getValue("seed:min_score");
00591
00592 for (Size s = min_spectra_; s < end_of_iteration; ++s)
00593 {
00594 ff_->setProgress(s);
00595
00596
00597 for (Size p = 0; p < map_[s].size(); ++p)
00598 {
00599 FloatDataArrays& meta = map_[s].getFloatDataArrays();
00600 DoubleReal overall_score = std::pow(meta[0][p] * meta[1][p] * meta[meta_index_isotope][p], 1.0f / 3.0f);
00601 meta[meta_index_overall][p] = overall_score;
00602
00603
00604 if (meta[2][p] != 0.0)
00605 {
00606
00607 if (!user_seeds && overall_score >= min_seed_score)
00608 {
00609 Seed seed;
00610 seed.spectrum = s;
00611 seed.peak = p;
00612 seed.intensity = map_[s][p].getIntensity();
00613 seeds.push_back(seed);
00614 }
00615
00616 else if (user_seeds && overall_score >= user_seed_score)
00617 {
00618
00619 FeatureType tmp;
00620 tmp.setMZ(map_[s][p].getMZ() - user_mz_tol);
00621 for (typename FeatureMapType::const_iterator it = std::lower_bound(seeds_.begin(), seeds_.end(), tmp, typename FeatureType::MZLess()); it < seeds_.end(); ++it)
00622 {
00623 if (it->getMZ() > map_[s][p].getMZ() + user_mz_tol)
00624 {
00625 break;
00626 }
00627 if (fabs(it->getMZ() - map_[s][p].getMZ()) < user_mz_tol &&
00628 fabs(it->getRT() - map_[s].getRT()) < user_rt_tol)
00629 {
00630 Seed seed;
00631 seed.spectrum = s;
00632 seed.peak = p;
00633 seed.intensity = map_[s][p].getIntensity();
00634 seeds.push_back(seed);
00635 break;
00636 }
00637 }
00638 }
00639 }
00640 }
00641 }
00642
00643 std::sort(seeds.rbegin(), seeds.rend());
00644
00645 if (debug_)
00646 {
00647
00648 FeatureMap<> seed_map;
00649 seed_map.reserve(seeds.size());
00650 for (Size i = 0; i < seeds.size(); ++i)
00651 {
00652 Size spectrum = seeds[i].spectrum;
00653 Size peak = seeds[i].peak;
00654 const FloatDataArrays& meta = map_[spectrum].getFloatDataArrays();
00655 Feature tmp;
00656 tmp.setIntensity(seeds[i].intensity);
00657 tmp.setOverallQuality(meta[meta_index_overall][peak]);
00658 tmp.setRT(map_[spectrum].getRT());
00659 tmp.setMZ(map_[spectrum][peak].getMZ());
00660 tmp.setMetaValue("intensity_score", meta[1][peak]);
00661 tmp.setMetaValue("pattern_score", meta[meta_index_isotope][peak]);
00662 tmp.setMetaValue("trace_score", meta[0][peak]);
00663 seed_map.push_back(tmp);
00664 }
00665 FeatureXMLFile().store(String("debug/seeds_") + String(c) + ".featureXML", seed_map);
00666 }
00667
00668 ff_->endProgress();
00669 std::cout << "Found " << seeds.size() << " seeds for charge " << c << "." << std::endl;
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683
00684 std::map<int, std::vector<int> > seeds_in_features;
00685 typedef std::map<int, OpenMS::Feature> FeatureMapType;
00686 FeatureMapType tmp_feature_map;
00687 gl_progress = 0;
00688 ff_->startProgress(0, seeds.size(), String("Extending seeds for charge ") + String(c));
00689 #ifdef _OPENMP
00690 #pragma omp parallel for
00691 #endif
00692 for (SignedSize i = 0; i < (SignedSize)seeds.size(); ++i)
00693 {
00694
00695
00696
00697
00698
00699 const SpectrumType& spectrum = map_[seeds[i].spectrum];
00700 const PeakType& peak = spectrum[seeds[i].peak];
00701
00702 IF_MASTERTHREAD
00703 {
00704 ff_->setProgress(gl_progress++);
00705
00706 if (debug_)
00707 {
00708 log_ << std::endl << "Seed " << i << ":" << std::endl;
00709
00710 log_ << " - Int: " << peak.getIntensity() << std::endl;
00711 log_ << " - RT: " << spectrum.getRT() << std::endl;
00712 log_ << " - MZ: " << peak.getMZ() << std::endl;
00713 }
00714 }
00715
00716
00717
00718 IsotopePattern best_pattern(0);
00719 DoubleReal isotope_fit_quality = findBestIsotopeFit_(seeds[i], c, best_pattern);
00720
00721 if (isotope_fit_quality < min_isotope_fit_)
00722 {
00723 abort_(seeds[i], "Could not find good enough isotope pattern containing the seed");
00724
00725 }
00726 else
00727 {
00728
00729
00730 MassTraces traces;
00731 traces.reserve(best_pattern.peak.size());
00732 extendMassTraces_(best_pattern, traces, meta_index_overall);
00733
00734
00735 DoubleReal seed_mz = map_[seeds[i].spectrum][seeds[i].peak].getMZ();
00736
00737 if (!traces.isValid(seed_mz, trace_tolerance_))
00738 {
00739 abort_(seeds[i], "Could not extend seed");
00740
00741 }
00742 else
00743 {
00744
00745
00746
00747
00748
00749 Int plot_nr = -1;
00750
00751 #ifdef _OPENMP
00752 #pragma omp critical (FeatureFinderAlgorithmPicked_PLOTNR)
00753 #endif
00754 {
00755 plot_nr = ++plot_nr_global;
00756 }
00757
00758
00759
00760
00761
00762 traces.updateBaseline();
00763 traces.baseline = 0.75 * traces.baseline;
00764
00765 traces[traces.max_trace].updateMaximum();
00766
00767
00768 double egh_tau = 0.0;
00769 TraceFitter<PeakType>* fitter = chooseTraceFitter_(egh_tau);
00770
00771 fitter->setParameters(trace_fitter_params);
00772 fitter->fit(traces);
00773
00774 #if 0
00775 TraceFitter<PeakType>* alt_fitter = new GaussTraceFitter<PeakType>();
00776 Param alt_p;
00777 alt_p.setValue("max_iteration", max_iterations);
00778 alt_p.setValue("epsilon_abs", epsilon_abs);
00779 alt_p.setValue("epsilon_rel", epsilon_rel);
00780
00781 alt_fitter->setParameters(alt_p);
00782 alt_fitter->fit(traces);
00783
00784 LOG_DEBUG << "EGH: " << fitter->getCenter() << " " << fitter->getHeight() << std::endl;
00785 LOG_DEBUG << "GAUSS: " << alt_fitter->getCenter() << " " << alt_fitter->getHeight() << std::endl;
00786 #endif
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799 MassTraces new_traces;
00800 cropFeature_(fitter, traces, new_traces);
00801
00802
00803
00804
00805
00806 String error_msg = "";
00807
00808 DoubleReal fit_score = 0.0;
00809 DoubleReal correlation = 0.0;
00810 DoubleReal final_score = 0.0;
00811
00812 bool feature_ok = checkFeatureQuality_(fitter, new_traces, seed_mz, min_feature_score, error_msg, fit_score, correlation, final_score);
00813 #ifdef _OPENMP
00814 #pragma omp critical (FeatureFinderAlgorithmPicked_DEBUG)
00815 #endif
00816 {
00817
00818 if (debug_)
00819 {
00820 writeFeatureDebugInfo_(fitter, traces, new_traces, feature_ok, error_msg, final_score, plot_nr, peak);
00821 }
00822 }
00823 traces = new_traces;
00824
00825
00826
00827 if (!feature_ok)
00828 {
00829 abort_(seeds[i], error_msg);
00830
00831 }
00832 else
00833 {
00834
00835
00836
00837
00838
00839 Feature f;
00840
00841 f.setMetaValue(3, plot_nr);
00842 f.setCharge(c);
00843 f.setOverallQuality(final_score);
00844 f.setMetaValue("score_fit", fit_score);
00845 f.setMetaValue("score_correlation", correlation);
00846 f.setRT(fitter->getCenter());
00847 f.setWidth(fitter->getFWHM());
00848
00849
00850 if (egh_tau != 0.0)
00851 {
00852 egh_tau = (static_cast<EGHTraceFitter<PeakType>*>(fitter))->getTau();
00853 f.setMetaValue("EGH_tau", egh_tau);
00854 f.setMetaValue("EGH_height", (static_cast<EGHTraceFitter<PeakType>*>(fitter))->getHeight());
00855 f.setMetaValue("EGH_sigma", (static_cast<EGHTraceFitter<PeakType>*>(fitter))->getSigmaSquare());
00856 }
00857
00858
00859 if (reported_mz_ == "maximum")
00860 {
00861 f.setMZ(traces[traces.getTheoreticalmaxPosition()].getAvgMZ());
00862 }
00863 else if (reported_mz_ == "average")
00864 {
00865 DoubleReal total_intensity = 0.0;
00866 DoubleReal average_mz = 0.0;
00867 for (Size t = 0; t < traces.size(); ++t)
00868 {
00869 for (Size p = 0; p < traces[t].peaks.size(); ++p)
00870 {
00871 average_mz += traces[t].peaks[p].second->getMZ() * traces[t].peaks[p].second->getIntensity();
00872 total_intensity += traces[t].peaks[p].second->getIntensity();
00873 }
00874 }
00875 average_mz /= total_intensity;
00876 f.setMZ(average_mz);
00877 }
00878 else if (reported_mz_ == "monoisotopic")
00879 {
00880 DoubleReal mono_mz = traces[traces.getTheoreticalmaxPosition()].getAvgMZ();
00881 mono_mz -= (Constants::PROTON_MASS_U / c) * (traces.getTheoreticalmaxPosition() + best_pattern.theoretical_pattern.trimmed_left);
00882 f.setMZ(mono_mz);
00883 }
00884
00885
00886
00887
00888 f.setIntensity(
00889 fitter->getFeatureIntensityContribution()
00890 / getIsotopeDistribution_(f.getMZ()).max);
00891
00892
00893 delete fitter;
00894
00895
00896 for (Size j = 0; j < traces.size(); ++j)
00897 {
00898 f.getConvexHulls().push_back(traces[j].getConvexhull());
00899 }
00900
00901 #ifdef _OPENMP
00902 #pragma omp critical (FeatureFinderAlgorithmPicked_TMPFEATUREMAP)
00903 #endif
00904 {
00905 tmp_feature_map[i] = f;
00906 }
00907
00908
00909
00910 DBoundingBox<2> bb = f.getConvexHull().getBoundingBox();
00911 for (Size j = i + 1; j < seeds.size(); ++j)
00912 {
00913 DoubleReal rt = map_[seeds[j].spectrum].getRT();
00914 DoubleReal mz = map_[seeds[j].spectrum][seeds[j].peak].getMZ();
00915 if (bb.encloses(rt, mz) && f.encloses(rt, mz))
00916 {
00917 #ifdef _OPENMP
00918 #pragma omp critical (FeatureFinderAlgorithmPicked_SEEDSINFEATURES)
00919 #endif
00920 {
00921 seeds_in_features[i].push_back(j);
00922 }
00923 }
00924 }
00925 }
00926 }
00927 }
00928 }
00929
00930
00931
00932
00933
00934 std::vector<int> seeds_contained;
00935 for (typename std::map<int, FeatureType>::iterator iter = tmp_feature_map.begin(); iter != tmp_feature_map.end(); ++iter)
00936 {
00937 int seed_nr = iter->first;
00938 bool is_used = false;
00939 for (Size i = 0; i < seeds_contained.size(); ++i)
00940 {
00941 if (seed_nr == seeds_contained[i]) { is_used = true; break; }
00942 }
00943 if (!is_used)
00944 {
00945 ++feature_candidates;
00946
00947
00948 iter->second.setMetaValue(3, feature_nr_global);
00949 ++feature_nr_global;
00950 features_->push_back(iter->second);
00951
00952 std::vector<int> curr_seed = seeds_in_features[seed_nr];
00953 for (Size k = 0; k < curr_seed.size(); ++k)
00954 {
00955 seeds_contained.push_back(curr_seed[k]);
00956
00957 }
00958 }
00959 }
00960
00961 IF_MASTERTHREAD ff_->endProgress();
00962 std::cout << "Found " << feature_candidates << " feature candidates for charge " << c << "." << std::endl;
00963 }
00964
00965
00966
00967
00968
00969
00970 ff_->startProgress(0, features_->size() * features_->size(), "Resolving overlapping features");
00971 if (debug_) log_ << "Resolving intersecting features (" << features_->size() << " candidates)" << std::endl;
00972
00973 features_->sortByMZ();
00974
00975 std::vector<DBoundingBox<2> > bbs(features_->size());
00976 DoubleReal max_mz_span = 0.0;
00977
00978 for (Size i = 0; i < features_->size(); ++i)
00979 {
00980 bbs[i] = (*features_)[i].getConvexHull().getBoundingBox();
00981 if (bbs[i].height() > max_mz_span)
00982 {
00983 max_mz_span = bbs[i].height();
00984 }
00985 }
00986
00987 Size removed(0);
00988
00989 for (Size i = 0; i < features_->size(); ++i)
00990 {
00991 Feature& f1((*features_)[i]);
00992 for (Size j = i + 1; j < features_->size(); ++j)
00993 {
00994 ff_->setProgress(i * features_->size() + j);
00995 Feature& f2((*features_)[j]);
00996
00997 if (f2.getMZ() - f1.getMZ() > 2.0 * max_mz_span) break;
00998
00999 if (f1.getIntensity() == 0.0 || f2.getIntensity() == 0.0) continue;
01000
01001 if (!bbs[i].intersects(bbs[j])) continue;
01002
01003 DoubleReal intersection = intersection_(f1, f2);
01004
01005 if (intersection >= max_feature_intersection_)
01006 {
01007 ++removed;
01008
01009 if (debug_) log_ << " - Intersection (" << (i + 1) << "/" << (j + 1) << "): " << intersection << std::endl;
01010 if (f1.getCharge() == f2.getCharge())
01011 {
01012 if (f1.getIntensity() * f1.getOverallQuality() > f2.getIntensity() * f2.getOverallQuality())
01013 {
01014 if (debug_) log_ << " - same charge -> removing duplicate " << (j + 1) << std::endl;
01015 f1.getSubordinates().push_back(f2);
01016 f2.setIntensity(0.0);
01017 }
01018 else
01019 {
01020 if (debug_) log_ << " - same charge -> removing duplicate " << (i + 1) << std::endl;
01021 f2.getSubordinates().push_back(f1);
01022 f1.setIntensity(0.0);
01023 }
01024 }
01025 else if (f2.getCharge() % f1.getCharge() == 0)
01026 {
01027 if (debug_) log_ << " - different charge (one is the multiple of the other) -> removing lower charge " << (i + 1) << std::endl;
01028 f2.getSubordinates().push_back(f1);
01029 f1.setIntensity(0.0);
01030 }
01031 else if (f1.getCharge() % f2.getCharge() == 0)
01032 {
01033 if (debug_) log_ << " - different charge (one is the multiple of the other) -> removing lower charge " << (i + 1) << std::endl;
01034 f1.getSubordinates().push_back(f2);
01035 f2.setIntensity(0.0);
01036 }
01037 else
01038 {
01039 if (f1.getOverallQuality() > f2.getOverallQuality())
01040 {
01041 if (debug_) log_ << " - different charge -> removing lower score " << (j + 1) << std::endl;
01042 f1.getSubordinates().push_back(f2);
01043 f2.setIntensity(0.0);
01044 }
01045 else
01046 {
01047 if (debug_) log_ << " - different charge -> removing lower score " << (i + 1) << std::endl;
01048 f2.getSubordinates().push_back(f1);
01049 f1.setIntensity(0.0);
01050 }
01051 }
01052 }
01053 }
01054 }
01055 LOG_INFO << "Removed " << removed << " overlapping features." << std::endl;
01056
01057 FeatureMap<> tmp;
01058 tmp.reserve(features_->size());
01059 for (Size i = 0; i < features_->size(); ++i)
01060 {
01061 if (features_->operator[](i).getIntensity() != 0.0)
01062 {
01063 tmp.push_back(features_->operator[](i));
01064 }
01065 }
01066 tmp.Base::swap(*features_);
01067
01068 features_->sortByIntensity(true);
01069 ff_->endProgress();
01070 std::cout << features_->size() << " features left." << std::endl;
01071
01072
01073 std::cout << std::endl;
01074 std::cout << "Abort reasons during feature construction:" << std::endl;
01075 for (std::map<String, UInt>::const_iterator it = aborts_.begin(); it != aborts_.end(); ++it)
01076 {
01077 std::cout << "- " << it->first << ": " << it->second << std::endl;
01078 }
01079 if (debug_)
01080 {
01081
01082 FeatureMap<> abort_map;
01083 abort_map.reserve(abort_reasons_.size());
01084 Size counter = 0;
01085 for (typename std::map<Seed, String>::iterator it2 = abort_reasons_.begin(); it2 != abort_reasons_.end(); ++it2, ++counter)
01086 {
01087 Feature f;
01088 f.setRT(map_[it2->first.spectrum].getRT());
01089 f.setMZ(map_[it2->first.spectrum][it2->first.peak].getMZ());
01090 f.setIntensity(map_[it2->first.spectrum][it2->first.peak].getIntensity());
01091 f.setMetaValue("label", it2->second);
01092 f.setUniqueId(counter);
01093 abort_map.push_back(f);
01094 }
01095 abort_map.setUniqueId();
01096 FeatureXMLFile().store("debug/abort_reasons.featureXML", abort_map);
01097
01098
01099 for (Size s = 0; s < map_.size(); ++s)
01100 {
01101 map_[s].getFloatDataArrays().erase(map_[s].getFloatDataArrays().begin() + 2);
01102 }
01103 MzMLFile().store("debug/input.mzML", map_);
01104 }
01105
01106 }
01107
01108 static FeatureFinderAlgorithm<PeakType, FeatureType>* create()
01109 {
01110 return new FeatureFinderAlgorithmPicked();
01111 }
01112
01113 static const String getProductName()
01114 {
01115 return "centroided";
01116 }
01117
01118 protected:
01120 MapType map_;
01122 mutable std::ofstream log_;
01124 bool debug_;
01126 std::map<String, UInt> aborts_;
01128 std::map<Seed, String> abort_reasons_;
01130 FeatureMapType seeds_;
01131
01133
01134 DoubleReal pattern_tolerance_;
01135 DoubleReal trace_tolerance_;
01136 UInt min_spectra_;
01137 UInt max_missing_trace_peaks_;
01138 DoubleReal slope_bound_;
01139 DoubleReal intensity_percentage_;
01140 DoubleReal intensity_percentage_optional_;
01141 DoubleReal optional_fit_improvement_;
01142 DoubleReal mass_window_width_;
01143 UInt intensity_bins_;
01144 DoubleReal min_isotope_fit_;
01145 DoubleReal min_trace_score_;
01146 DoubleReal min_rt_span_;
01147 DoubleReal max_rt_span_;
01148 DoubleReal max_feature_intersection_;
01149 String reported_mz_;
01150
01151
01153
01154
01155 DoubleReal intensity_rt_step_;
01157 DoubleReal intensity_mz_step_;
01159 std::vector<std::vector<std::vector<DoubleReal> > > intensity_thresholds_;
01161
01163 std::vector<TheoreticalIsotopePattern> isotope_distributions_;
01164
01165
01166 virtual void updateMembers_()
01167 {
01168 pattern_tolerance_ = param_.getValue("mass_trace:mz_tolerance");
01169 trace_tolerance_ = param_.getValue("isotopic_pattern:mz_tolerance");
01170 min_spectra_ = (UInt) std::floor((DoubleReal)param_.getValue("mass_trace:min_spectra") * 0.5);
01171 max_missing_trace_peaks_ = param_.getValue("mass_trace:max_missing");
01172 slope_bound_ = param_.getValue("mass_trace:slope_bound");
01173 intensity_percentage_ = (DoubleReal)param_.getValue("isotopic_pattern:intensity_percentage") / 100.0;
01174 intensity_percentage_optional_ = (DoubleReal)param_.getValue("isotopic_pattern:intensity_percentage_optional") / 100.0;
01175 optional_fit_improvement_ = (DoubleReal)param_.getValue("isotopic_pattern:optional_fit_improvement") / 100.0;
01176 mass_window_width_ = param_.getValue("isotopic_pattern:mass_window_width");
01177 intensity_bins_ = param_.getValue("intensity:bins");
01178 min_isotope_fit_ = param_.getValue("feature:min_isotope_fit");
01179 min_trace_score_ = param_.getValue("feature:min_trace_score");
01180 min_rt_span_ = param_.getValue("feature:min_rt_span");
01181 max_rt_span_ = param_.getValue("feature:max_rt_span");
01182 max_feature_intersection_ = param_.getValue("feature:max_intersection");
01183 reported_mz_ = param_.getValue("feature:reported_mz");
01184 }
01185
01187 void abort_(const Seed& seed, const String& reason)
01188 {
01189 if (debug_) log_ << "Abort: " << reason << std::endl;
01190 aborts_[reason]++;
01191 if (debug_) abort_reasons_[seed] = reason;
01192 }
01193
01198 DoubleReal intersection_(const Feature& f1, const Feature& f2) const
01199 {
01200
01201 DoubleReal s1 = 0.0;
01202 const std::vector<ConvexHull2D>& hulls1 = f1.getConvexHulls();
01203 for (Size i = 0; i < hulls1.size(); ++i)
01204 {
01205 s1 += hulls1[i].getBoundingBox().width();
01206 }
01207
01208
01209 DoubleReal s2 = 0.0;
01210 const std::vector<ConvexHull2D>& hulls2 = f2.getConvexHulls();
01211 for (Size j = 0; j < hulls2.size(); ++j)
01212 {
01213 s2 += hulls2[j].getBoundingBox().width();
01214 }
01215
01216
01217 DoubleReal overlap = 0.0;
01218 for (Size i = 0; i < hulls1.size(); ++i)
01219 {
01220 DBoundingBox<2> bb1 = hulls1[i].getBoundingBox();
01221 for (Size j = 0; j < hulls2.size(); ++j)
01222 {
01223 DBoundingBox<2> bb2 = hulls2[j].getBoundingBox();
01224 if (bb1.intersects(bb2))
01225 {
01226 if (bb1.minPosition()[0] <= bb2.minPosition()[0] &&
01227 bb1.maxPosition()[0] >= bb2.maxPosition()[0])
01228 {
01229 overlap += bb2.width();
01230 }
01231 else if (bb2.minPosition()[0] <= bb1.minPosition()[0] &&
01232 bb2.maxPosition()[0] >= bb1.maxPosition()[0])
01233 {
01234 overlap += bb1.width();
01235 }
01236 else if (bb1.minPosition()[0] <= bb2.minPosition()[0] &&
01237 bb1.maxPosition()[0] <= bb2.maxPosition()[0])
01238 {
01239 overlap += bb1.maxPosition()[0] - bb2.minPosition()[0];
01240 }
01241 else if (bb2.minPosition()[0] <= bb1.minPosition()[0] &&
01242 bb2.maxPosition()[0] <= bb1.maxPosition()[0])
01243 {
01244 overlap += bb2.maxPosition()[0] - bb1.minPosition()[0];
01245 }
01246 }
01247 }
01248 }
01249
01250 return overlap / std::min(s1, s2);
01251 }
01252
01254 const TheoreticalIsotopePattern& getIsotopeDistribution_(DoubleReal mass) const
01255 {
01256
01257 Size index = (Size) std::floor(mass / mass_window_width_);
01258
01259 if (index >= isotope_distributions_.size())
01260 {
01261 throw Exception::InvalidValue(__FILE__, __LINE__, __PRETTY_FUNCTION__, "IsotopeDistribution not precalculated. Maximum allowed index is " + String(isotope_distributions_.size()), String(index));
01262 }
01263
01264
01265 return isotope_distributions_[index];
01266 }
01267
01275 DoubleReal findBestIsotopeFit_(const Seed& center, UInt charge, IsotopePattern& best_pattern) const
01276 {
01277 if (debug_) log_ << "Testing isotope patterns for charge " << charge << ": " << std::endl;
01278 const SpectrumType& spectrum = map_[center.spectrum];
01279 const TheoreticalIsotopePattern& isotopes = getIsotopeDistribution_(spectrum[center.peak].getMZ() * charge);
01280 if (debug_) log_ << " - Seed: " << center.peak << " (mz:" << spectrum[center.peak].getMZ() << ")" << std::endl;
01281
01282
01283 DoubleReal mass_window = (DoubleReal)(isotopes.size() + 1) / (DoubleReal)charge;
01284 if (debug_) log_ << " - Mass window: " << mass_window << std::endl;
01285 Size end = center.peak;
01286 while (end < spectrum.size() &&
01287 spectrum[end].getMZ() < spectrum[center.peak].getMZ() + mass_window)
01288 {
01289 ++end;
01290 }
01291 --end;
01292
01293
01294 SignedSize begin = center.peak;
01295 while (begin >= 0 &&
01296 spectrum[begin].getMZ() > spectrum[center.peak].getMZ() - mass_window)
01297 {
01298 --begin;
01299 }
01300 ++begin;
01301 if (debug_) log_ << " - Begin: " << begin << " (mz:" << spectrum[begin].getMZ() << ")" << std::endl;
01302 if (debug_) log_ << " - End: " << end << " (mz:" << spectrum[end].getMZ() << ")" << std::endl;
01303
01304
01305 DoubleReal max_score = 0.0;
01306 for (Size start = begin; start <= end; ++start)
01307 {
01308
01309 Size peak_index = start;
01310 IsotopePattern pattern(isotopes.size());
01311 if (debug_) log_ << " - Fitting at " << start << " (mz:" << spectrum[start].getMZ() << ")" << std::endl;
01312 for (Size iso = 0; iso < isotopes.size(); ++iso)
01313 {
01314 DoubleReal pos = spectrum[start].getMZ() + iso / (DoubleReal)charge;
01315 findIsotope_(pos, center.spectrum, pattern, iso, peak_index);
01316 }
01317
01318
01319 bool seed_contained = false;
01320 for (Size iso = 0; iso < pattern.peak.size(); ++iso)
01321 {
01322 if (pattern.peak[iso] == (Int)center.peak && pattern.spectrum[iso] == center.spectrum)
01323 {
01324 seed_contained = true;
01325 break;
01326 }
01327 }
01328 if (!seed_contained)
01329 {
01330 if (debug_) log_ << " - aborting: seed is not contained!" << std::endl;
01331 continue;
01332 }
01333
01334 DoubleReal score = isotopeScore_(isotopes, pattern, false);
01335
01336
01337 seed_contained = false;
01338 for (Size iso = 0; iso < pattern.peak.size(); ++iso)
01339 {
01340 if (pattern.peak[iso] == (Int)center.peak &&
01341 pattern.spectrum[iso] == center.spectrum)
01342 {
01343 seed_contained = true;
01344 break;
01345 }
01346 }
01347 if (!seed_contained)
01348 {
01349 if (debug_) log_ << " - aborting: seed was removed during isotope fit!" << std::endl;
01350 continue;
01351 }
01352
01353 if (debug_) log_ << " - final score: " << score << std::endl;
01354 if (score > max_score)
01355 {
01356 max_score = score;
01357 best_pattern = pattern;
01358 }
01359 }
01360 if (debug_) log_ << " - best score : " << max_score << std::endl;
01361 best_pattern.theoretical_pattern = isotopes;
01362 return max_score;
01363 }
01364
01372 void extendMassTraces_(const IsotopePattern& pattern, MassTraces& traces, Size meta_index_overall) const
01373 {
01374
01375 DoubleReal max_int = 0.0;
01376 Size max_trace_index = 0;
01377 for (Size p = 0; p < pattern.peak.size(); ++p)
01378 {
01379 if (pattern.peak[p] < 0) continue;
01380 if (map_[pattern.spectrum[p]][pattern.peak[p]].getIntensity() > max_int)
01381 {
01382 max_int = map_[pattern.spectrum[p]][pattern.peak[p]].getIntensity();
01383 max_trace_index = p;
01384 }
01385 }
01386
01387
01388 Size start_index = pattern.spectrum[max_trace_index];
01389 const PeakType* start_peak = &(map_[pattern.spectrum[max_trace_index]][pattern.peak[max_trace_index]]);
01390 DoubleReal start_mz = start_peak->getMZ();
01391 DoubleReal start_rt = map_[start_index].getRT();
01392 if (debug_) log_ << " - Trace " << max_trace_index << " (maximum intensity)" << std::endl;
01393 if (debug_) log_ << " - extending from: " << map_[start_index].getRT() << " / " << start_mz << " (int: " << start_peak->getIntensity() << ")" << std::endl;
01394
01395 MassTrace max_trace;
01396 max_trace.peaks.push_back(std::make_pair(start_rt, start_peak));
01397 extendMassTrace_(max_trace, start_index, start_mz, false, meta_index_overall);
01398 extendMassTrace_(max_trace, start_index, start_mz, true, meta_index_overall);
01399
01400 DoubleReal rt_max = max_trace.peaks.back().first;
01401 DoubleReal rt_min = max_trace.peaks.begin()->first;
01402 if (debug_) log_ << " - rt bounds: " << rt_min << "-" << rt_max << std::endl;
01403
01404 if (!max_trace.isValid() || max_trace.peaks.size() < 2 * min_spectra_ - max_missing_trace_peaks_)
01405 {
01406 if (debug_) log_ << " - could not extend trace with maximum intensity => abort" << std::endl;
01407 return;
01408 }
01409 for (Size p = 0; p < pattern.peak.size(); ++p)
01410 {
01411 if (debug_) log_ << " - Trace " << p << std::endl;
01412 if (p == max_trace_index)
01413 {
01414 if (debug_) log_ << " - previously extended maximum trace" << std::endl;
01415 traces.push_back(max_trace);
01416 traces.back().theoretical_int = pattern.theoretical_pattern.intensity[p];
01417 traces.max_trace = traces.size() - 1;
01418 continue;
01419 }
01420 Seed starting_peak;
01421 starting_peak.spectrum = pattern.spectrum[p];
01422 starting_peak.peak = pattern.peak[p];
01423 if (pattern.peak[p] == -2)
01424 {
01425 if (debug_) log_ << " - removed during isotope fit" << std::endl;
01426 continue;
01427 }
01428 else if (pattern.peak[p] == -1)
01429 {
01430 if (debug_) log_ << " - missing" << std::endl;
01431 continue;
01432 }
01433 starting_peak.intensity = map_[starting_peak.spectrum][starting_peak.peak].getIntensity();
01434 if (debug_) log_ << " - trace seed: " << map_[starting_peak.spectrum].getRT() << " / " << map_[starting_peak.spectrum][starting_peak.peak].getMZ() << " (int: " << map_[starting_peak.spectrum][starting_peak.peak].getIntensity() << ")" << std::endl;
01435
01436
01437 Size begin = std::max((Size)0, starting_peak.spectrum - min_spectra_);
01438 Size end = std::min(starting_peak.spectrum + min_spectra_, (Size)map_.size());
01439 DoubleReal mz = map_[starting_peak.spectrum][starting_peak.peak].getMZ();
01440 DoubleReal inte = map_[starting_peak.spectrum][starting_peak.peak].getIntensity();
01441 for (Size spectrum_index = begin; spectrum_index < end; ++spectrum_index)
01442 {
01443
01444 SignedSize peak_index = -1;
01445 try
01446 {
01447 peak_index = map_[spectrum_index].findNearest(map_[starting_peak.spectrum][starting_peak.peak].getMZ());
01448 }
01449 catch (...)
01450 {
01451 peak_index = -1;
01452 }
01453
01454 if (peak_index < 0 ||
01455 map_[spectrum_index][peak_index].getIntensity() <= inte ||
01456 std::fabs(mz - map_[spectrum_index][peak_index].getMZ()) >= pattern_tolerance_
01457 )
01458 {
01459 continue;
01460 }
01461
01462 starting_peak.spectrum = spectrum_index;
01463 starting_peak.peak = peak_index;
01464 inte = map_[spectrum_index][peak_index].getIntensity();
01465 }
01466 if (debug_) log_ << " - extending from: " << map_[starting_peak.spectrum].getRT() << " / " << map_[starting_peak.spectrum][starting_peak.peak].getMZ() << " (int: " << map_[starting_peak.spectrum][starting_peak.peak].getIntensity() << ")" << std::endl;
01467
01468
01469
01470 MassTrace trace;
01471 const PeakType* seed = &(map_[starting_peak.spectrum][starting_peak.peak]);
01472
01473 trace.peaks.push_back(std::make_pair(map_[starting_peak.spectrum].getRT(), seed));
01474 extendMassTrace_(trace, starting_peak.spectrum, seed->getMZ(), false, meta_index_overall, rt_min, rt_max);
01475 extendMassTrace_(trace, starting_peak.spectrum, seed->getMZ(), true, meta_index_overall, rt_min, rt_max);
01476
01477
01478 if (!trace.isValid())
01479 {
01480 if (debug_) log_ << " - could not extend trace " << std::endl;
01481
01482 if (p < traces.max_trace)
01483 {
01484 traces.clear();
01485 continue;
01486 }
01487 else if (p > traces.max_trace)
01488 {
01489 break;
01490 }
01491 }
01492 traces.push_back(trace);
01493 traces.back().theoretical_int = pattern.theoretical_pattern.intensity[p];
01494 }
01495 }
01496
01515 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
01516 {
01517
01518 if (increase_rt)
01519 {
01520 ++spectrum_index;
01521 std::reverse(trace.peaks.begin(), trace.peaks.end());
01522 }
01523 else
01524 {
01525 --spectrum_index;
01526 }
01527
01528
01529 bool boundaries = false;
01530 if (max_rt != min_rt)
01531 {
01532 boundaries = true;
01533 }
01534
01535
01536 DoubleReal current_slope_bound = (1.0 + (DoubleReal)boundaries) * slope_bound_;
01537
01538 Size delta_count = min_spectra_;
01539 std::vector<DoubleReal> deltas(delta_count - 1, 0);
01540
01541 DoubleReal last_observed_intensity = trace.peaks.back().second->getIntensity();
01542
01543 UInt missing_peaks = 0;
01544 Size peaks_before_extension = trace.peaks.size();
01545 String abort_reason = "";
01546
01547 while ((!increase_rt && spectrum_index >= 0) || (increase_rt && spectrum_index < (SignedSize)map_.size()))
01548 {
01549 if (boundaries &&
01550 ((!increase_rt && map_[spectrum_index].getRT() < min_rt) ||
01551 (increase_rt && map_[spectrum_index].getRT() > max_rt))
01552 )
01553 {
01554 abort_reason = "Hit upper/lower boundary";
01555 break;
01556 }
01557
01558 SignedSize peak_index = -1;
01559
01560 try
01561 {
01562 peak_index = map_[spectrum_index].findNearest(mz);
01563 }
01564 catch (...)
01565 {
01566 peak_index = -1;
01567 }
01568
01569
01570 if (
01571 peak_index < 0
01572 || map_[spectrum_index].getFloatDataArrays()[meta_index_overall][peak_index] < 0.01
01573 || positionScore_(mz, map_[spectrum_index][peak_index].getMZ(), trace_tolerance_) == 0.0
01574 )
01575 {
01576 ++missing_peaks;
01577
01578 if (missing_peaks > max_missing_trace_peaks_)
01579 {
01580 abort_reason = "too many peaks missing";
01581 break;
01582 }
01583 }
01584 else
01585 {
01586 missing_peaks = 0;
01587
01588
01589 trace.peaks.push_back(std::make_pair(map_[spectrum_index].getRT(), &(map_[spectrum_index][peak_index])));
01590
01591
01592 deltas.push_back((map_[spectrum_index][peak_index].getIntensity() - last_observed_intensity) / last_observed_intensity);
01593 last_observed_intensity = map_[spectrum_index][peak_index].getIntensity();
01594
01595
01596 DoubleReal average_delta = std::accumulate(deltas.end() - delta_count, deltas.end(), 0.0) / (DoubleReal)delta_count;
01597 if (average_delta > current_slope_bound)
01598 {
01599 abort_reason = String("Average delta above threshold: ") + average_delta + "/" + current_slope_bound;
01600
01601
01602 Size remove = std::min((Size)(trace.peaks.size() - peaks_before_extension), delta_count - 1);
01603 trace.peaks.erase(trace.peaks.end() - remove, trace.peaks.end());
01604 break;
01605 }
01606 }
01607
01608
01609 if (increase_rt) ++spectrum_index;
01610 else --spectrum_index;
01611 }
01612 if (debug_) log_ << " - Added " << (trace.peaks.size() - peaks_before_extension) << " peaks (abort: " << abort_reason << ")" << std::endl;
01613 }
01614
01616 template <typename SpectrumType>
01617 Size nearest_(DoubleReal pos, const SpectrumType& spec, Size start) const
01618 {
01619 Size index = start;
01620 DoubleReal distance = std::fabs(pos - spec[index].getMZ());
01621 ++index;
01622 while (index < spec.size())
01623 {
01624 DoubleReal new_distance = std::fabs(pos - spec[index].getMZ());
01625 if (new_distance < distance)
01626 {
01627 distance = new_distance;
01628 ++index;
01629 }
01630 else
01631 {
01632 break;
01633 }
01634 }
01635 return --index;
01636 }
01637
01647 void findIsotope_(DoubleReal pos, Size spectrum_index, IsotopePattern& pattern, Size pattern_index, Size& peak_index) const
01648 {
01649 if (debug_) log_ << " - Isotope " << pattern_index << ": ";
01650
01651 DoubleReal intensity = 0.0;
01652 DoubleReal pos_score = 0.0;
01653 UInt matches = 0;
01654
01655
01656 const SpectrumType& spectrum = map_[spectrum_index];
01657 peak_index = nearest_(pos, spectrum, peak_index);
01658 DoubleReal mz_score = positionScore_(pos, spectrum[peak_index].getMZ(), pattern_tolerance_);
01659 pattern.theoretical_mz[pattern_index] = pos;
01660
01661 if (mz_score != 0.0)
01662 {
01663 if (debug_) log_ << String::number(spectrum[peak_index].getIntensity(), 1) << " ";
01664 pattern.peak[pattern_index] = peak_index;
01665 pattern.spectrum[pattern_index] = spectrum_index;
01666 intensity += spectrum[peak_index].getIntensity();
01667 pos_score += mz_score;
01668 ++matches;
01669 }
01670
01671
01672 if (spectrum_index != 0 && !map_[spectrum_index - 1].empty())
01673 {
01674 const SpectrumType& spectrum_before = map_[spectrum_index - 1];
01675 Size index_before = spectrum_before.findNearest(pos);
01676 DoubleReal mz_score = positionScore_(pos, spectrum_before[index_before].getMZ(), pattern_tolerance_);
01677 if (mz_score != 0.0)
01678 {
01679 if (debug_) log_ << String::number(spectrum_before[index_before].getIntensity(), 1) << "b ";
01680 intensity += spectrum_before[index_before].getIntensity();
01681 pos_score += mz_score;
01682 ++matches;
01683
01684 if (pattern.peak[pattern_index] == -1)
01685 {
01686 pattern.peak[pattern_index] = index_before;
01687 pattern.spectrum[pattern_index] = spectrum_index - 1;
01688 }
01689 }
01690 }
01691
01692 if (spectrum_index != map_.size() - 1 && !map_[spectrum_index + 1].empty())
01693 {
01694 const SpectrumType& spectrum_after = map_[spectrum_index + 1];
01695 Size index_after = spectrum_after.findNearest(pos);
01696 DoubleReal mz_score = positionScore_(pos, spectrum_after[index_after].getMZ(), pattern_tolerance_);
01697 if (mz_score != 0.0)
01698 {
01699 if (debug_) log_ << String::number(spectrum_after[index_after].getIntensity(), 1) << "a ";
01700 intensity += spectrum_after[index_after].getIntensity();
01701 pos_score += mz_score;
01702 ++matches;
01703
01704 if (pattern.peak[pattern_index] == -1)
01705 {
01706 pattern.peak[pattern_index] = index_after;
01707 pattern.spectrum[pattern_index] = spectrum_index + 1;
01708 }
01709 }
01710 }
01711
01712 if (matches == 0)
01713 {
01714 if (debug_) log_ << " missing" << std::endl;
01715 pattern.peak[pattern_index] = -1;
01716 pattern.mz_score[pattern_index] = 0.0;
01717 pattern.intensity[pattern_index] = 0.0;
01718 }
01719 else
01720 {
01721 if (debug_) log_ << "=> " << intensity / matches << std::endl;
01722 pattern.mz_score[pattern_index] = pos_score / matches;
01723 pattern.intensity[pattern_index] = intensity / matches;
01724 }
01725 }
01726
01728 DoubleReal positionScore_(DoubleReal pos1, DoubleReal pos2, DoubleReal allowed_deviation) const
01729 {
01730 DoubleReal diff = fabs(pos1 - pos2);
01731 if (diff <= 0.5 * allowed_deviation)
01732 {
01733 return 0.1 * (0.5 * allowed_deviation - diff) / (0.5 * allowed_deviation) + 0.9;
01734 }
01735 else if (diff <= allowed_deviation)
01736 {
01737 return 0.9 * (allowed_deviation - diff) / (0.5 * allowed_deviation);
01738 }
01739 return 0.0;
01740 }
01741
01743 DoubleReal isotopeScore_(const TheoreticalIsotopePattern& isotopes, IsotopePattern& pattern, bool consider_mz_distances) const
01744 {
01745 if (debug_) log_ << " - fitting " << pattern.intensity.size() << " peaks" << std::endl;
01746
01747 for (Size iso = 0 + isotopes.optional_begin; iso < pattern.peak.size() - isotopes.optional_end; ++iso)
01748 {
01749 if (pattern.peak[iso] == -1)
01750 {
01751 if (debug_) log_ << " - aborting: core peak is missing" << std::endl;
01752 return 0.0;
01753 }
01754 }
01755
01756
01757
01758 DoubleReal best_int_score = 0.01;
01759 Size best_begin = 0;
01760 for (Size i = isotopes.optional_begin; i > 0; --i)
01761 {
01762 if (pattern.peak[i - 1] == -1)
01763 {
01764 best_begin = i;
01765 break;
01766 }
01767 }
01768 Size best_end = 0;
01769 for (Size i = isotopes.optional_end; i > 0; --i)
01770 {
01771 if (pattern.peak[pattern.peak.size() - i] == -1)
01772 {
01773 best_end = i;
01774 break;
01775 }
01776 }
01777 if (debug_) log_ << " - best_begin/end: " << best_begin << "/" << best_end << std::endl;
01778 for (Size b = best_begin; b <= isotopes.optional_begin; ++b)
01779 {
01780 for (Size e = best_end; e <= isotopes.optional_end; ++e)
01781 {
01782
01783 if (isotopes.size() - b - e > 2 || (b == best_begin &&
01784 e == best_end &&
01785 isotopes.size() - b - e > 1))
01786 {
01787 DoubleReal int_score = Math::pearsonCorrelationCoefficient(isotopes.intensity.begin() + b, isotopes.intensity.end() - e, pattern.intensity.begin() + b, pattern.intensity.end() - e);
01788 if (boost::math::isnan(int_score)) int_score = 0.0;
01789 if (isotopes.size() - b - e == 2 && int_score > min_isotope_fit_) int_score = min_isotope_fit_;
01790 if (debug_) log_ << " - fit (" << b << "/" << e << "): " << int_score;
01791 if (int_score / best_int_score >= 1.0 + optional_fit_improvement_)
01792 {
01793 if (debug_) log_ << " - new best fit ";
01794 best_int_score = int_score;
01795 best_begin = b;
01796 best_end = e;
01797 }
01798 if (debug_) log_ << std::endl;
01799 }
01800 }
01801 }
01802
01803
01804 if (pattern.mz_score.size() - best_begin - best_end == 0)
01805 {
01806 return 0.0;
01807 }
01808
01809
01810 for (Size i = 0; i < best_begin; ++i)
01811 {
01812 pattern.peak[i] = -2;
01813 pattern.intensity[i] = 0.0;
01814 pattern.mz_score[i] = 0.0;
01815 }
01816
01817 for (Size i = 0; i < best_end; ++i)
01818 {
01819 pattern.peak[isotopes.size() - 1 - i] = -2;
01820 pattern.intensity[isotopes.size() - 1 - i] = 0.0;
01821 pattern.mz_score[isotopes.size() - 1 - i] = 0.0;
01822 }
01823
01824 if (consider_mz_distances)
01825 {
01826 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);
01827 }
01828
01829
01830 OPENMS_POSTCONDITION(best_int_score >= 0.0, (String("Internal error: Isotope score (") + best_int_score + ") should be >=0.0").c_str())
01831 OPENMS_POSTCONDITION(best_int_score <= 1.0, (String("Internal error: Isotope score (") + best_int_score + ") should be <=1.0").c_str())
01832 return best_int_score;
01833 }
01834
01845 DoubleReal intensityScore_(Size spectrum, Size peak) const
01846 {
01847
01848 DoubleReal intensity = map_[spectrum][peak].getIntensity();
01849 DoubleReal rt = map_[spectrum].getRT();
01850 DoubleReal mz = map_[spectrum][peak].getMZ();
01851 DoubleReal rt_min = map_.getMinRT();
01852 DoubleReal mz_min = map_.getMinMZ();
01853 UInt rt_bin = std::min(2 * intensity_bins_ - 1, (UInt) std::floor((rt - rt_min) / intensity_rt_step_ * 2.0));
01854 UInt mz_bin = std::min(2 * intensity_bins_ - 1, (UInt) std::floor((mz - mz_min) / intensity_mz_step_ * 2.0));
01855
01856 UInt ml, mh;
01857 if (mz_bin == 0 || mz_bin == 2 * intensity_bins_ - 1)
01858 {
01859 ml = mz_bin / 2;
01860 mh = mz_bin / 2;
01861 }
01862 else if (Math::isOdd(mz_bin))
01863 {
01864 ml = mz_bin / 2;
01865 mh = mz_bin / 2 + 1;
01866 }
01867 else
01868 {
01869 ml = mz_bin / 2 - 1;
01870 mh = mz_bin / 2;
01871 }
01872
01873 UInt rl, rh;
01874 if (rt_bin == 0 || rt_bin == 2 * intensity_bins_ - 1)
01875 {
01876 rl = rt_bin / 2;
01877 rh = rt_bin / 2;
01878 }
01879 else if (Math::isOdd(rt_bin))
01880 {
01881 rl = rt_bin / 2;
01882 rh = rt_bin / 2 + 1;
01883 }
01884 else
01885 {
01886 rl = rt_bin / 2 - 1;
01887 rh = rt_bin / 2;
01888 }
01889
01890 DoubleReal drl = std::fabs(rt_min + (0.5 + rl) * intensity_rt_step_ - rt) / intensity_rt_step_;
01891 DoubleReal drh = std::fabs(rt_min + (0.5 + rh) * intensity_rt_step_ - rt) / intensity_rt_step_;
01892 DoubleReal dml = std::fabs(mz_min + (0.5 + ml) * intensity_mz_step_ - mz) / intensity_mz_step_;
01893 DoubleReal dmh = std::fabs(mz_min + (0.5 + mh) * intensity_mz_step_ - mz) / intensity_mz_step_;
01894
01895
01896 DoubleReal d1 = std::sqrt(std::pow(1.0 - drl, 2) + std::pow(1.0 - dml, 2));
01897 DoubleReal d2 = std::sqrt(std::pow(1.0 - drh, 2) + std::pow(1.0 - dml, 2));
01898 DoubleReal d3 = std::sqrt(std::pow(1.0 - drl, 2) + std::pow(1.0 - dmh, 2));
01899 DoubleReal d4 = std::sqrt(std::pow(1.0 - drh, 2) + std::pow(1.0 - dmh, 2));
01900 DoubleReal d_sum = d1 + d2 + d3 + d4;
01901
01902
01903 DoubleReal final = intensityScore_(rl, ml, intensity) * (d1 / d_sum)
01904 + intensityScore_(rh, ml, intensity) * (d2 / d_sum)
01905 + intensityScore_(rl, mh, intensity) * (d3 / d_sum)
01906 + intensityScore_(rh, mh, intensity) * (d4 / d_sum);
01907
01908 OPENMS_POSTCONDITION(final >= 0.0, (String("Internal error: Intensity score (") + final + ") should be >=0.0").c_str())
01909 OPENMS_POSTCONDITION(final <= 1.0001, (String("Internal error: Intensity score (") + final + ") should be <=1.0").c_str())
01910 return final;
01911 }
01912
01919 TraceFitter<PeakType>* chooseTraceFitter_(double& tau)
01920 {
01921
01922 if (param_.getValue("feature:rt_shape") == "asymmetric")
01923 {
01924 LOG_DEBUG << "use asymmetric rt peak shape" << std::endl;
01925 tau = -1.0;
01926 return new EGHTraceFitter<PeakType>();
01927 }
01928 else
01929 {
01930 LOG_DEBUG << "use symmetric rt peak shape" << std::endl;
01931 return new GaussTraceFitter<PeakType>();
01932 }
01933 }
01934
01935 DoubleReal intensityScore_(Size rt_bin, Size mz_bin, DoubleReal intensity) const
01936 {
01937
01938 const std::vector<DoubleReal>& quantiles20 = intensity_thresholds_[rt_bin][mz_bin];
01939
01940 std::vector<DoubleReal>::const_iterator it = std::lower_bound(quantiles20.begin(), quantiles20.end(), intensity);
01941
01942 if (it == quantiles20.end())
01943 {
01944 return 1.0;
01945 }
01946
01947 DoubleReal bin_score = 0.0;
01948 if (it == quantiles20.begin())
01949 {
01950 bin_score = 0.05 * intensity / *it;
01951 }
01952 else
01953 {
01954
01955 bin_score = 0.05 * (intensity - *(it - 1)) / (*it - *(it - 1));
01956 }
01957
01958 DoubleReal final = bin_score +
01959 0.05 * ((it - quantiles20.begin()) - 1.0);
01960
01961
01962 if (final < 0.0) final = 0.0;
01963 if (final > 1.0) final = 1.0;
01964
01965
01966 return final;
01967 }
01968
01975
01984 void cropFeature_(TraceFitter<PeakType>* fitter,
01985 const MassTraces& traces,
01986 MassTraces& new_traces)
01987 {
01988 DoubleReal low_bound = fitter->getLowerRTBound();
01989 DoubleReal high_bound = fitter->getUpperRTBound();
01990
01991 if (debug_) log_ << " => RT bounds: " << low_bound << " - " << high_bound << std::endl;
01992 for (Size t = 0; t < traces.size(); ++t)
01993 {
01994 const MassTrace& trace = traces[t];
01995 if (debug_) log_ << " - Trace " << t << ": (" << trace.theoretical_int << ")" << std::endl;
01996
01997 MassTrace new_trace;
01998
01999 DoubleReal deviation = 0.0;
02000 std::vector<DoubleReal> v_theo, v_real;
02001 for (Size k = 0; k < trace.peaks.size(); ++k)
02002 {
02003
02004 if (trace.peaks[k].first >= low_bound && trace.peaks[k].first <= high_bound)
02005 {
02006 new_trace.peaks.push_back(trace.peaks[k]);
02007
02008 DoubleReal theo = traces.baseline + fitter->computeTheoretical(trace, k);
02009
02010 v_theo.push_back(theo);
02011 DoubleReal real = trace.peaks[k].second->getIntensity();
02012 v_real.push_back(real);
02013 deviation += std::fabs(real - theo) / theo;
02014 }
02015 }
02016 DoubleReal fit_score = 0.0;
02017 DoubleReal correlation = 0.0;
02018 DoubleReal final_score = 0.0;
02019 if (!new_trace.peaks.empty())
02020 {
02021 fit_score = deviation / new_trace.peaks.size();
02022 correlation = std::max(0.0, Math::pearsonCorrelationCoefficient(v_theo.begin(), v_theo.end(), v_real.begin(), v_real.end()));
02023 final_score = std::sqrt(correlation * std::max(0.0, 1.0 - fit_score));
02024 }
02025 if (debug_) log_ << " - peaks: " << new_trace.peaks.size() << " / " << trace.peaks.size() << " - relative deviation: " << fit_score << " - correlation: " << correlation << " - final score: " << correlation << std::endl;
02026
02027 if (!new_trace.isValid() || final_score < min_trace_score_)
02028 {
02029 if (t < traces.max_trace)
02030 {
02031 new_traces = MassTraces();
02032 if (debug_) log_ << " - removed this and previous traces due to bad fit" << std::endl;
02033 new_traces.clear();
02034 continue;
02035 }
02036 else if (t == traces.max_trace)
02037 {
02038 new_traces = MassTraces();
02039 if (debug_) log_ << " - aborting (max trace was removed)" << std::endl;
02040 break;
02041 }
02042 else if (t > traces.max_trace)
02043 {
02044 if (debug_) log_ << " - removed due to bad fit => omitting the rest" << std::endl;
02045 break;
02046 }
02047 }
02048
02049 else
02050 {
02051 new_trace.theoretical_int = trace.theoretical_int;
02052 new_traces.push_back(new_trace);
02053 if (t == traces.max_trace)
02054 {
02055 new_traces.max_trace = new_traces.size() - 1;
02056 }
02057 }
02058 }
02059 new_traces.baseline = traces.baseline;
02060 }
02061
02085 bool checkFeatureQuality_(TraceFitter<PeakType>* fitter,
02086 MassTraces& feature_traces,
02087 const DoubleReal& seed_mz, const DoubleReal& min_feature_score,
02088 String& error_msg, DoubleReal& fit_score, DoubleReal& correlation, DoubleReal& final_score)
02089 {
02090 bool feature_ok = true;
02091
02092
02093 if (feature_ok)
02094 {
02095
02096 if (fitter->checkMaximalRTSpan(max_rt_span_))
02097 {
02098 feature_ok = false;
02099 error_msg = "Invalid fit: Fitted model is bigger than 'max_rt_span'";
02100 }
02101 }
02102
02103
02104 if (!feature_traces.isValid(seed_mz, trace_tolerance_))
02105 {
02106 feature_ok = false;
02107 error_msg = "Invalid feature after fit - too few traces or peaks left";
02108 }
02109
02110
02111 if (feature_ok)
02112 {
02113 std::pair<DoubleReal, DoubleReal> rt_bounds = feature_traces.getRTBounds();
02114 if (fitter->getCenter() < rt_bounds.first || fitter->getCenter() > rt_bounds.second)
02115 {
02116 feature_ok = false;
02117 error_msg = "Invalid fit: Center outside of feature bounds";
02118 }
02119 }
02120
02121
02122 if (feature_ok)
02123 {
02124 std::pair<DoubleReal, DoubleReal> rt_bounds = feature_traces.getRTBounds();
02125 if (fitter->checkMinimalRTSpan(rt_bounds, min_rt_span_))
02126 {
02127 feature_ok = false;
02128 error_msg = "Invalid fit: Less than 'min_rt_span' left after fit";
02129 }
02130 }
02131
02132
02133 if (feature_ok)
02134 {
02135 std::vector<DoubleReal> v_theo, v_real;
02136 DoubleReal deviation = 0.0;
02137 for (Size t = 0; t < feature_traces.size(); ++t)
02138 {
02139 MassTrace& trace = feature_traces[t];
02140 for (Size k = 0; k < trace.peaks.size(); ++k)
02141 {
02142
02143 DoubleReal theo = feature_traces.baseline + fitter->computeTheoretical(trace, k);
02144 v_theo.push_back(theo);
02145 DoubleReal real = trace.peaks[k].second->getIntensity();
02146 v_real.push_back(real);
02147 deviation += std::fabs(real - theo) / theo;
02148 }
02149 }
02150 fit_score = std::max(0.0, 1.0 - (deviation / feature_traces.getPeakCount()));
02151 correlation = std::max(0.0, Math::pearsonCorrelationCoefficient(v_theo.begin(), v_theo.end(), v_real.begin(), v_real.end()));
02152 final_score = std::sqrt(correlation * fit_score);
02153
02154 if (final_score < min_feature_score)
02155 {
02156 feature_ok = false;
02157 error_msg = "Feature quality too low after fit";
02158 }
02159
02160
02161 if (debug_)
02162 {
02163 log_ << "Quality estimation:" << std::endl;
02164 log_ << " - relative deviation: " << fit_score << std::endl;
02165 log_ << " - correlation: " << correlation << std::endl;
02166 log_ << " => final score: " << final_score << std::endl;
02167 }
02168 }
02169
02170 return feature_ok;
02171 }
02172
02186 void writeFeatureDebugInfo_(TraceFitter<PeakType>* fitter,
02187 const MassTraces& traces,
02188 const MassTraces& new_traces,
02189 bool feature_ok, const String error_msg, const DoubleReal final_score, const Int plot_nr, const PeakType& peak,
02190 const String path = "debug/features/")
02191 {
02192
02193 DoubleReal pseudo_rt_shift = param_.getValue("debug:pseudo_rt_shift");
02194 TextFile tf;
02195
02196 String script = String("plot \"") + path + plot_nr + ".dta\" title 'before fit (RT: " + String::number(fitter->getCenter(), 2) + " m/z: " + String::number(peak.getMZ(), 4) + ")' with points 1";
02197
02198 for (Size k = 0; k < traces.size(); ++k)
02199 {
02200 for (Size j = 0; j < traces[k].peaks.size(); ++j)
02201 {
02202 tf.push_back(String(pseudo_rt_shift * k + traces[k].peaks[j].first) + "\t"+ traces[k].peaks[j].second->getIntensity());
02203 }
02204 }
02205 tf.store(path + plot_nr + ".dta");
02206
02207 if (new_traces.getPeakCount() != 0)
02208 {
02209 tf.clear();
02210 for (Size k = 0; k < new_traces.size(); ++k)
02211 {
02212 for (Size j = 0; j < new_traces[k].peaks.size(); ++j)
02213 {
02214 tf.push_back(String(pseudo_rt_shift * k + new_traces[k].peaks[j].first) + "\t" + new_traces[k].peaks[j].second->getIntensity());
02215 }
02216 }
02217
02218 tf.store(path + plot_nr + "_cropped.dta");
02219 script = script + ", \"" + path + plot_nr + "_cropped.dta\" title 'feature ";
02220
02221 if (!feature_ok)
02222 {
02223 script = script + " - " + error_msg;
02224 }
02225 else
02226 {
02227 script = script + (features_->size() + 1) + " (score: " + String::number(final_score, 3) + ")";
02228 }
02229 script = script + "' with points 3";
02230 }
02231
02232 tf.clear();
02233 for (Size k = 0; k < traces.size(); ++k)
02234 {
02235 char fun = 'f';
02236 fun += (char)k;
02237 tf.push_back(fitter->getGnuplotFormula(traces[k], fun, traces.baseline, pseudo_rt_shift * k));
02238
02239 script = script + ", " + fun + "(x) title 'Trace " + k + " (m/z: " + String::number(traces[k].getAvgMZ(), 4) + ")'";
02240 }
02241
02242
02243 tf.push_back("set xlabel \"pseudo RT (mass traces side-by-side)\"");
02244 tf.push_back("set ylabel \"intensity\"");
02245 tf.push_back("set samples 1000");
02246 tf.push_back(script);
02247 tf.push_back("pause -1");
02248 tf.store(path + plot_nr + ".plot");
02249 }
02250
02252 private:
02253
02255 FeatureFinderAlgorithmPicked& operator=(const FeatureFinderAlgorithmPicked&);
02257 FeatureFinderAlgorithmPicked(const FeatureFinderAlgorithmPicked&);
02258 };
02259
02260 }
02261
02262 #endif // OPENMS_TRANSFORMATIONS_FEATUREFINDER_FEATUREFINDERALGORITHMPICKED_H