36 #ifndef OPENMS_FILTERING_CALIBRATION_INTERNALCALIBRATION_H
37 #define OPENMS_FILTERING_CALIBRATION_INTERNALCALIBRATION_H
46 #include <gsl/gsl_fit.h>
85 template <
typename InputPeakType>
100 template <
typename InputPeakType>
116 template <
typename InputPeakType>
138 void calibrateMapGlobally(
const FeatureMap<> & feature_map,
FeatureMap<> & calibrated_feature_map, std::vector<PeptideIdentification> & ref_ids,
String trafo_file_name =
"");
142 template <
typename InputPeakType>
149 void makeLinearRegression_(std::vector<DoubleReal> & observed_masses, std::vector<DoubleReal> & theoretical_masses);
152 void checkReferenceIds_(std::vector<PeptideIdentification> & pep_ids);
155 void checkReferenceIds_(
const FeatureMap<> & feature_map);
165 template <
typename InputPeakType>
168 #ifdef DEBUG_CALIBRATION
169 std::cout.precision(writtenDigits<DoubleReal>());
173 std::cout <<
"Input is empty." << std::endl;
179 std::cout <<
"Attention: this function is assuming peak data." << std::endl;
181 calibrated_exp = exp;
183 Size num_ref_peaks = ref_masses.size();
184 bool use_ppm = (
param_.
getValue(
"mz_tolerance_unit") ==
"ppm") ?
true :
false;
188 for (
Size spec = 0; spec < exp.
size(); ++spec)
191 if (exp[spec].getMSLevel() != 1)
197 std::vector<DoubleReal> corr_masses, rel_errors, found_ref_masses;
199 for (
Size peak = 0; peak < exp[spec].
size(); ++peak)
201 for (
Size ref_peak = 0; ref_peak < num_ref_peaks; ++ref_peak)
203 if (!use_ppm && fabs(exp[spec][peak].getMZ() - ref_masses[ref_peak]) < mz_tol)
205 found_ref_masses.push_back(ref_masses[ref_peak]);
206 corr_masses.push_back(exp[spec][peak].getMZ());
210 else if (use_ppm && fabs(exp[spec][peak].getMZ() - ref_masses[ref_peak]) / ref_masses[ref_peak] * 1e6 < mz_tol)
212 found_ref_masses.push_back(ref_masses[ref_peak]);
213 corr_masses.push_back(exp[spec][peak].getMZ());
221 std::cout <<
"spec: " << spec
222 <<
" less than 2 reference masses were detected within a reasonable error range\n";
223 std::cout <<
"This spectrum cannot be calibrated!\n";
228 for (
Size ref_peak = 0; ref_peak < found_ref_masses.size(); ++ref_peak)
230 rel_errors.push_back((found_ref_masses[ref_peak] - corr_masses[ref_peak]) / corr_masses[ref_peak] * 1e6);
236 for (
unsigned int peak = 0; peak < calibrated_exp[spec].
size(); ++peak)
238 #ifdef DEBUG_CALIBRATION
239 std::cout << calibrated_exp[spec][peak].getMZ() <<
"\t";
241 DoubleReal mz = calibrated_exp[spec][peak].getMZ();
243 calibrated_exp[spec][peak].setMZ(mz);
244 #ifdef DEBUG_CALIBRATION
245 std::cout << calibrated_exp[spec][peak].getMZ() << std::endl;
254 template <
typename InputPeakType>
257 std::vector<PeptideIdentification> & ref_ids,
String trafo_file_name)
259 bool use_ppm =
param_.
getValue(
"mz_tolerance_unit") ==
"ppm" ?
true :
false;
263 std::cout <<
"Input is empty." << std::endl;
269 std::cout <<
"Attention: this function is assuming peak data." << std::endl;
274 std::vector<DoubleReal> theoretical_masses, observed_masses;
275 for (
Size p_id = 0; p_id < ref_ids.size(); ++p_id)
277 for (
Size p_h = 0; p_h < ref_ids[p_id].getHits().size(); ++p_h)
279 Int charge = ref_ids[p_id].getHits()[p_h].getCharge();
283 while (rt_iter != exp.
begin() && rt_iter->getMSLevel() != 1)
292 if ((mz_iter + 1) != rt_iter->
end()
293 && fabs((mz_iter + 1)->getMZ() - (
DoubleReal)ref_ids[p_id].getMetaValue(
"MZ")) < fabs(dist)
294 && mz_iter != rt_iter->
begin()
295 && fabs((mz_iter - 1)->getMZ() - (
DoubleReal)ref_ids[p_id].getMetaValue(
"MZ")) < fabs((mz_iter + 1)->getMZ() - (
DoubleReal)ref_ids[p_id].getMetaValue(
"MZ")))
298 fabs((mz_iter + 1)->getMZ() - (
DoubleReal)ref_ids[p_id].getMetaValue(
"MZ")) / (
DoubleReal)ref_ids[p_id].getMetaValue(
"MZ") * 1e06 < mz_tolerance) ||
299 (!use_ppm && fabs((mz_iter + 1)->getMZ() - (
DoubleReal)ref_ids[p_id].getMetaValue(
"MZ")) < mz_tolerance))
302 observed_masses.push_back((mz_iter + 1)->getMZ());
303 theoretical_masses.push_back(theo_mass);
308 else if (mz_iter != rt_iter->
begin()
309 && fabs((mz_iter - 1)->getMZ() - (
DoubleReal)ref_ids[p_id].getMetaValue(
"MZ")) < fabs(dist))
312 fabs((mz_iter - 1)->getMZ() - (
DoubleReal)ref_ids[p_id].getMetaValue(
"MZ")) / (
DoubleReal)ref_ids[p_id].getMetaValue(
"MZ") * 1e06 < mz_tolerance) ||
313 (!use_ppm && fabs((mz_iter - 1)->getMZ() - (
DoubleReal)ref_ids[p_id].getMetaValue(
"MZ")) < mz_tolerance))
316 observed_masses.push_back((mz_iter - 1)->getMZ());
317 theoretical_masses.push_back(theo_mass);
325 fabs((mz_iter)->getMZ() - (
DoubleReal)ref_ids[p_id].getMetaValue(
"MZ")) / (
DoubleReal)ref_ids[p_id].getMetaValue(
"MZ") * 1e06 < mz_tolerance) ||
326 (!use_ppm && fabs((mz_iter)->getMZ() - (
DoubleReal)ref_ids[p_id].getMetaValue(
"MZ")) < mz_tolerance))
329 observed_masses.push_back(mz_iter->getMZ());
330 theoretical_masses.push_back(theo_mass);
343 for (
Size spec = 0; spec < exp.
size(); ++spec)
346 if (exp[spec].getMSLevel() != 1)
348 calibrated_exp[spec] = exp[spec];
352 calibrated_exp[spec] = exp[spec];
354 for (
unsigned int peak = 0; peak < exp[spec].
size(); ++peak)
356 #ifdef DEBUG_CALIBRATION
357 std::cout << exp[spec][peak].getMZ() <<
"\t";
361 calibrated_exp[spec][peak].setMZ(mz);
362 #ifdef DEBUG_CALIBRATION
363 std::cout << calibrated_exp[spec][peak].getMZ() << std::endl;
368 if (trafo_file_name !=
"")
374 template <
typename InputPeakType>
379 std::cout <<
"Input is empty." << std::endl;
385 std::cout <<
"Attention: this function is assuming peak data." << std::endl;
389 Size num_ref_peaks = ref_masses.size();
390 bool use_ppm = (
param_.
getValue(
"mz_tolerance_unit") ==
"ppm") ?
true :
false;
393 std::vector<DoubleReal> corr_masses, rel_errors, found_ref_masses;
396 for (
Size spec = 0; spec < exp.
size(); ++spec)
399 if (exp[spec].getMSLevel() != 1)
401 for (
Size peak = 0; peak < exp[spec].
size(); ++peak)
403 for (
Size ref_peak = 0; ref_peak < num_ref_peaks; ++ref_peak)
405 if (!use_ppm && fabs(exp[spec][peak].getMZ() - ref_masses[ref_peak]) < mz_tol)
407 found_ref_masses.push_back(ref_masses[ref_peak]);
408 corr_masses.push_back(exp[spec][peak].getMZ());
412 else if (use_ppm && fabs(exp[spec][peak].getMZ() - ref_masses[ref_peak]) / ref_masses[ref_peak] * 1e6 < mz_tol)
414 found_ref_masses.push_back(ref_masses[ref_peak]);
415 corr_masses.push_back(exp[spec][peak].getMZ());
424 std::cout <<
"Less than 2 reference masses were detected within a reasonable error range\n";
425 std::cout <<
"This spectrum cannot be calibrated!\n";
435 for (
Size spec = 0; spec < exp.
size(); ++spec)
438 if (exp[spec].getMSLevel() != 1)
440 calibrated_exp[spec] = exp[spec];
445 calibrated_exp[spec] = exp[spec];
447 for (
unsigned int peak = 0; peak < exp[spec].
size(); ++peak)
449 #ifdef DEBUG_CALIBRATION
450 std::cout << exp[spec][peak].getMZ() <<
"\t";
454 calibrated_exp[spec][peak].setMZ(mz);
456 #ifdef DEBUG_CALIBRATION
457 std::cout << calibrated_exp[spec][peak].getMZ() << std::endl;
464 if (trafo_file_name !=
"")
472 #endif // OPENMS_FILTERING_CALIBRATION_INTERNALCALIBRATION_H
void calibrateMapGlobally(const MSExperiment< InputPeakType > &exp, MSExperiment< InputPeakType > &calibrated_exp, std::vector< DoubleReal > &ref_masses, String trafo_file_name="")
Calibrate a peak map using given reference masses with one calibration function for the whole map...
Definition: InternalCalibration.h:375
A more convenient string class.
Definition: String.h:56
Size size() const
Definition: MSExperiment.h:117
Peak data (also called centroided data or stick data)
Definition: SpectrumSettings.h:74
Param param_
Container for current parameters.
Definition: DefaultParamHandler.h:148
void makeLinearRegression_(std::vector< DoubleReal > &observed_masses, std::vector< DoubleReal > &theoretical_masses)
the actual calibration function
A simple calibration method using linear interpolation of given reference masses. ...
Definition: InternalCalibration.h:62
ContainerType::const_iterator ConstIterator
Non-mutable iterator.
Definition: MSSpectrum.h:125
Iterator begin()
Definition: MSExperiment.h:147
void resize(Size s)
Definition: MSExperiment.h:122
Iterator end()
Definition: MSExperiment.h:157
ConstIterator RTBegin(CoordinateType rt) const
Fast search for spectrum range begin.
Definition: MSExperiment.h:349
void endProgress() const
Ends the progress display.
~InternalCalibration()
Destructor.
Definition: InternalCalibration.h:71
const DataValue & getValue(const String &key) const
Returns a value of a parameter.
Definition: Residue.h:363
Representation of a mass spectrometry experiment.
Definition: MSExperiment.h:68
void calibrateMapSpectrumwise(const MSExperiment< InputPeakType > &exp, MSExperiment< InputPeakType > &calibrated_exp, std::vector< DoubleReal > &ref_masses)
Calibrate a peak map using given reference masses with a separate calibration function for each spect...
Definition: InternalCalibration.h:166
std::vector< SpectrumType >::const_iterator ConstIterator
Non-mutable iterator.
Definition: MSExperiment.h:103
size_t Size
Size type e.g. used as variable which can hold result of size()
Definition: Types.h:144
void checkReferenceIds_(std::vector< PeptideIdentification > &pep_ids)
check if reference ids contain RT and MZ information as meta values
Base class for all classes that want to report their progess.
Definition: ProgressLogger.h:56
void startProgress(SignedSize begin, SignedSize end, const String &label) const
Initializes the progress display.
bool empty() const
Definition: MSExperiment.h:127
void setProgress(SignedSize value) const
Sets the current progress.
A base class for all classes handling default parameters.
Definition: DefaultParamHandler.h:90
int Int
Signed integer type.
Definition: Types.h:100
Description of the experimental settings.
Definition: ExperimentalSettings.h:59
TransformationDescription trafo_
here the transformation is stored
Definition: InternalCalibration.h:161
double DoubleReal
Double-precision real type.
Definition: Types.h:118