35 #ifndef OPENMS_TRANSFORMATIONS_FEATUREFINDER_FEATUREFINDERALGORITHMMRM_H
36 #define OPENMS_TRANSFORMATIONS_FEATUREFINDER_FEATUREFINDERALGORITHMMRM_H
49 #include <boost/math/special_functions/fpclassify.hpp>
72 template <
class PeakType,
class FeatureType>
104 defaults_.
setValue(
"min_rt_distance", 10.0,
"Minimal distance of MRM features in seconds.");
108 defaults_.
setValue(
"min_signal_to_noise_ratio", 2.0,
"Minimal S/N ratio a peak must have to be taken into account. Set to zero if the MRM-traces contains mostly signals, and no noise.");
113 defaults_.
setValue(
"resample_traces",
"false",
"If set to true, each trace, which is in this case a part of the MRM monitoring trace with signal is resampled, using the minimal distance of two data points in RT dimension",
StringList::create(
"advanced"));
149 std::cerr <<
"Starting feature finding #chromatograms=" <<
map_->
getChromatograms().size() <<
", #spectra=" <<
map_->
size() << std::endl;
152 typename std::vector<MSChromatogram<ChromatogramPeak> >::const_iterator first_it =
map_->
getChromatograms().begin();
161 peak.
setMZ(it->getRT());
163 chromatogram.push_back(peak);
176 DoubleReal min_distance(std::numeric_limits<DoubleReal>::max()), old_rt(0);
181 std::cerr <<
"CHROMATOGRAM: " << it->getMZ() <<
" " << it->getIntensity() << std::endl;
184 if (rt_diff < min_distance && rt_diff > 0)
186 min_distance = rt_diff;
188 old_rt = it->getMZ();
193 std::cerr <<
"Min_distance=" << min_distance << std::endl;
195 if (min_distance > 50 || chromatogram.size() < min_num_peaks_per_feature)
200 resampler_param.
setValue(
"spacing", min_distance);
202 resampler.
raster(chromatogram);
216 std::cerr <<
"win_len m/z: " << (chromatogram.end() - 1)->getMZ() <<
" " << chromatogram.begin()->getMZ() << std::endl;
218 sne_param.setValue(
"win_len", (chromatogram.end() - 1)->getMZ() - chromatogram.begin()->getMZ());
220 if ((
DoubleReal)sne_param.getValue(
"win_len") < 10e-4)
225 sne.
init(chromatogram.begin(), chromatogram.end());
229 std::cerr << first_it->getPrecursor().getMZ() <<
" " << first_it->getProduct().getMZ() <<
" ";
234 sit->setMetaValue(
"SN", sn);
237 std::cerr << sit->getMZ() <<
" " << sit->getIntensity() <<
" " << sn << std::endl;
239 if (min_signal_to_noise_ratio == 0 || sn > min_signal_to_noise_ratio)
241 sn_chrom.push_back(*sit);
247 std::vector<std::vector<DPosition<2> > > sections;
252 std::cerr <<
"SECTIONS: " << sit->getMZ() <<
" " << sit->getIntensity() << std::endl;
254 this_rt = sit->getMZ();
255 if (sections.empty() || (this_rt - last_rt) > min_rt_distance)
259 std::cerr <<
"Starting new section, sections.size()=" << sections.size() <<
", rt_diff=" << this_rt - last_rt << std::endl;
262 std::vector<DPosition<2> > section;
263 section.push_back(
DPosition<2>(this_rt, sit->getIntensity()));
264 sections.push_back(section);
268 sections.back().push_back(
DPosition<2>(this_rt, sit->getIntensity()));
275 for (
Size i = 0; i != sections.size(); ++i)
277 if (sections[i].size() < min_num_peaks_per_feature)
282 std::vector<DoubleReal> deltas(2, 0.0);
284 for (
Size j = 1; j < sections[i].size(); ++j)
286 deltas.push_back((sections[i][j].getY() - last_int) / last_int);
287 last_int = sections[i][j].getY();
290 DoubleReal average_delta = std::accumulate(deltas.begin(), deltas.end(), 0.0) / (
DoubleReal)3.0;
291 std::cerr <<
"AverageDelta: " << average_delta <<
" (" << sections[i][j].getX() <<
", " << sections[i][j].getY() <<
")" << std::endl;
298 for (
Size i = 0; i != sections.size(); ++i)
300 if (sections[i].size() > min_num_peaks_per_feature)
304 for (
Size j = 0; j != sections[i].size(); ++j)
307 p.
setMZ(sections[i][j].getX());
309 filter_spec.push_back(p);
314 std::vector<DoubleReal> distances;
315 for (
Size j = 1; j < filter_spec.size(); ++j)
317 distances.push_back(filter_spec[j].getMZ() - filter_spec[j - 1].getMZ());
319 DoubleReal dist_average = std::accumulate(distances.begin(), distances.end(), 0.0) / (
DoubleReal)distances.size();
324 new_peak.
setMZ(filter_spec.back().getMZ() + dist_average);
325 filter_spec.push_back(new_peak);
326 new_peak.
setMZ(filter_spec.back().getMZ() + dist_average);
327 filter_spec.push_back(new_peak);
328 new_peak.
setMZ(filter_spec.back().getMZ() + dist_average);
329 filter_spec.push_back(new_peak);
332 new_peak.
setMZ(filter_spec.front().getMZ() - dist_average);
333 filter_spec.insert(filter_spec.begin(), new_peak);
334 new_peak.
setMZ(filter_spec.front().getMZ() - dist_average);
335 filter_spec.insert(filter_spec.begin(), new_peak);
336 new_peak.
setMZ(filter_spec.front().getMZ() - dist_average);
337 filter_spec.insert(filter_spec.begin(), new_peak);
339 filter_param.setValue(
"gaussian_width", 4 * dist_average);
343 filter.
filter(filter_spec);
346 std::vector<PeakType> data_to_fit;
347 for (
Size j = 0; j != filter_spec.size(); ++j)
352 data_to_fit.push_back(p);
363 for (
Size j = 0; j < sections[i].size(); ++j)
365 hull_points[j][0] = sections[i][j].getX();
366 hull_points[j][1] = first_it->getProduct().getMZ();
368 rt_sum += sections[i][j].getX();
369 intensity_sum += sections[i][j].getY();
391 if (write_debug_files)
394 std::ofstream data_out(
String(base_name +
"_data.dat").c_str());
395 for (
Size j = 0; j < sections[i].size(); ++j)
400 data_out << rt <<
" " << intensity << std::endl;
404 std::ofstream smoothed_data_out(
String(base_name +
"_smoothed_data.dat").c_str());
405 for (
Size j = 0; j < filter_spec.size(); ++j)
407 smoothed_data_out << filter_spec[j].getMZ() + 0.5 <<
" " << filter_spec[j].getIntensity() << std::endl;
409 smoothed_data_out.close();
411 std::ofstream fit_out(
String(base_name +
"_rt_fit.dat").c_str());
418 for (
DoubleReal pos = bb_min; pos < bb_max; pos += int_step)
421 fit_out << pos <<
" " << emg_model.
getIntensity(pos) << std::endl;
425 std::ofstream gnuplot_out(
String(base_name +
"_gnuplot.gpl").c_str());
426 gnuplot_out <<
"set terminal png" << std::endl;
427 gnuplot_out <<
"set output \"" << base_name <<
".png\"" << std::endl;
428 gnuplot_out <<
"plot '" << base_name <<
"_data.dat' w i, '" << base_name <<
"_smoothed_data.dat' w i, '" << base_name <<
"_rt_fit.dat' w lp title 'quality=" << f.
getOverallQuality() <<
"'" << std::endl;
430 String gnuplot_call =
"gnuplot " + base_name +
"_gnuplot.gpl";
431 int error = system(gnuplot_call.c_str());
434 std::cerr <<
"An error occurred during the gnuplot execution" << std::endl;
478 quality = fitter.
fit1d(rt_input_data, model);
481 if (boost::math::isnan(quality)) quality = -1.0;
495 #endif // OPENMS_TRANSFORMATIONS_FEATUREFINDER_FEATUREFINDERALGORITHMMRM_H
virtual double getSignalToNoise(const PeakIterator &data_point)
Definition: SignalToNoiseEstimator.h:130
ProductModel & setModel(UInt dim, BaseModel< 1 > *dist)
set model dist for dimension dim
Definition: ProductModel.h:184
Exponentially modified gaussian distribution model for elution profiles.
Definition: EmgModel.h:50
FeatureMapType * features_
Output data pointer.
Definition: FeatureFinderAlgorithm.h:144
Param defaults_
Container for default parameters. This member should be filled in the constructor of derived classes!...
Definition: DefaultParamHandler.h:155
void setValue(const String &key, const DataValue &value, const String &description="", const StringList &tags=StringList())
Sets a value.
A more convenient string class.
Definition: String.h:56
Linear Resampling of raw data.
Definition: LinearResampler.h:61
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
FeatureFinderAlgorithm for MRM experiments.
Definition: FeatureFinderAlgorithmMRM.h:73
Size size() const
Definition: MSExperiment.h:117
void setMZ(CoordinateType coordinate)
Mutable access to the m/z coordinate (index 1)
Definition: Peak2D.h:197
QualityType getOverallQuality() const
Non-mutable access to the overall quality.
Param param_
Container for current parameters.
Definition: DefaultParamHandler.h:148
The representation of a chromatogram.
Definition: MSChromatogram.h:53
const MapType * map_
Input data pointer.
Definition: FeatureFinderAlgorithm.h:141
FeatureFinderAlgorithm< PeakType, FeatureType >::MapType MapType
Definition: FeatureFinderAlgorithmMRM.h:80
Retention time dimension id (0 if used as a const int)
Definition: Peak2D.h:76
std::vector< PointType > PointArrayType
Definition: ConvexHull2D.h:79
Abstract base class for FeatureFinder algorithms.
Definition: FeatureFinderAlgorithm.h:74
ContainerType::const_iterator ConstIterator
Non-mutable iterator.
Definition: MSSpectrum.h:125
MapType::SpectrumType SpectrumType
Definition: FeatureFinderAlgorithmMRM.h:81
Exponentially modified gaussian distribution fitter (1-dim.) using Levenberg-Marquardt algorithm (GSL...
Definition: EmgFitter1D.h:47
Mass-to-charge dimension id (1 if used as a const int)
Definition: Peak2D.h:77
QualityType fit1d(const RawDataArrayType &range, InterpolationModel *&model)
return interpolation model
void setModelDescription(const ModelDescription< 2 > &q)
Set the model description.
CoordinateType getRT() const
Returns the RT coordinate (index 0)
Definition: Peak2D.h:203
The class template is only implemented for D=2 because we use Peak2D here.
Definition: ProductModel.h:65
void setMZ(CoordinateType mz)
Mutable access to m/z.
Definition: Peak1D.h:114
void setIntensity(IntensityType intensity)
Mutable access to the data point intensity (height)
Definition: Peak1D.h:105
A 1-dimensional raw data point or peak mith meta information.
Definition: RichPeak1D.h:52
ContainerType::iterator Iterator
Mutable iterator.
Definition: MSSpectrum.h:123
void setParameters(const Param ¶m)
Sets the parameters.
A 2-dimensional hull representation in [counter]clockwise direction - depending on axis labelling...
Definition: ConvexHull2D.h:75
The purpose of this struct is to provide definitions of classes and typedefs which are used throughou...
Definition: FeatureFinderDefs.h:51
DoubleReal fitRT_(std::vector< PeakType > &rt_input_data, InterpolationModel *&model) const
Definition: FeatureFinderAlgorithmMRM.h:458
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
Definition: FeatureFinderAlgorithmMRM.h:95
Abstract class for 1D-models that are approximated using linear interpolation.
Definition: InterpolationModel.h:55
Estimates the signal/noise (S/N) ratio of each data point in a scan based on an iterative scheme whic...
Definition: SignalToNoiseEstimatorMeanIterative.h:69
virtual void init(const PeakIterator &it_begin, const PeakIterator &it_end)
Set the start and endpoint of the raw data interval, for which signal to noise ratios will be estimat...
Definition: SignalToNoiseEstimator.h:111
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
IntensityType getIntensity(const PositionType &pos) const
access model predicted intensity at position pos
Definition: InterpolationModel.h:102
const std::vector< ConvexHull2D > & getConvexHulls() const
Non-mutable access to the convex hulls.
void setQuality(Size index, QualityType q)
Set the quality in dimension c.
SpectrumType::FloatDataArrays FloatDataArrays
Definition: FeatureFinderAlgorithmMRM.h:82
A 1-dimensional raw data point or peak.
Definition: Peak1D.h:55
FeatureFinderAlgorithmMRM()
default constructor
Definition: FeatureFinderAlgorithmMRM.h:101
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 setPosition(const PositionType &position)
Mutable access to the position.
Definition: Peak2D.h:185
virtual void updateMembers_()
This method is used to update extra member variables at the end of the setParameters() method...
Definition: FeatureFinderAlgorithmMRM.h:487
An LC-MS feature.
Definition: Feature.h:66
virtual void run()
Main method for actual FeatureFinder.
Definition: FeatureFinderAlgorithmMRM.h:123
void setValidStrings(const String &key, const std::vector< String > &strings)
Sets the valid strings for the parameter key.
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
Representation of a mass spectrometry experiment.
Definition: MSExperiment.h:68
This class represents a Gaussian lowpass-filter which works on uniform as well as on non-uniform prof...
Definition: GaussFilter.h:73
void setOverallQuality(QualityType q)
Set the overall quality.
void setMinInt(const String &key, Int min)
Sets the minimum value for the integer or integer list parameter key.
void raster(MSSpectrum< PeakType > &spectrum)
Applies the resampling algorithm to an MSSpectrum.
Definition: LinearResampler.h:85
size_t Size
Size type e.g. used as variable which can hold result of size()
Definition: Types.h:144
void startProgress(SignedSize begin, SignedSize end, const String &label) const
Initializes the progress display.
bool toBool() const
Conversion to bool.
void setProgress(SignedSize value) const
Sets the current progress.
void addPoints(const PointArrayType &points)
Definition: FeatureFinderAlgorithmMRM.h:96
static const String getProductName()
Definition: FeatureFinderAlgorithmMRM.h:451
Map class based on the STL map (containing serveral convenience functions)
Definition: Map.h:50
static FeatureFinderAlgorithm< PeakType, FeatureType > * create()
Definition: FeatureFinderAlgorithmMRM.h:446
void setMinFloat(const String &key, DoubleReal min)
Sets the minimum value for the floating point or floating point list parameter key.
void filter(MSSpectrum< PeakType > &spectrum)
Smoothes an MSSpectrum containing profile data.
Definition: GaussFilter.h:92
const Param & getParameters() const
Non-mutable access to the parameters.
void setSamples()
set sample/supporting points of interpolation
const std::vector< MSChromatogram< ChromatogramPeakType > > & getChromatograms() const
returns the chromatogram list
Definition: MSExperiment.h:768
void defaultsToParam_()
Updates the parameters after the defaults have been set in the constructor.
double DoubleReal
Double-precision real type.
Definition: Types.h:118