35 #ifndef OPENMS_ANALYSIS_OPENSWATH_MRMTRANSITIONGROUPPICKER_H
36 #define OPENMS_ANALYSIS_OPENSWATH_MRMTRANSITIONGROUPPICKER_H
114 template <
typename SpectrumT,
typename TransitionT>
117 std::vector<RichPeakChromatogram> picked_chroms_;
131 RichPeakChromatogram picked_chrom;
134 picked_chroms_.push_back(picked_chrom);
141 int chr_idx, peak_idx, cnt = 0;
142 std::vector<MRMFeature> features;
145 chr_idx = -1; peak_idx = -1;
146 findLargestPeak(picked_chroms_, chr_idx, peak_idx);
147 if (chr_idx == -1 && peak_idx == -1)
break;
150 MRMFeature mrm_feature = createMRMFeature(transition_group, picked_chroms_, chr_idx, peak_idx);
153 features.push_back(mrm_feature);
157 if ((stop_after_feature_ > 0 && cnt > stop_after_feature_) &&
165 for (
Size i = 0; i < features.size(); i++)
169 for (
Size j = 0; j < i; j++)
171 if ((
double)mrm_feature.
getMetaValue(
"leftWidth") >= (
double)features[j].getMetaValue(
"leftWidth") &&
184 template <
typename SpectrumT,
typename TransitionT>
186 std::vector<SpectrumT>& picked_chroms,
int& chr_idx,
int& peak_idx)
190 double best_left = picked_chroms[chr_idx].getFloatDataArrays()[1][peak_idx];
191 double best_right = picked_chroms[chr_idx].getFloatDataArrays()[2][peak_idx];
192 double peak_apex = picked_chroms[chr_idx][peak_idx].getRT();
193 LOG_DEBUG <<
"**** Creating MRMFeature for peak " << chr_idx <<
" " << peak_idx <<
" " <<
194 picked_chroms[chr_idx][peak_idx] <<
" with borders " << best_left <<
" " <<
195 best_right <<
" (" << best_right - best_left <<
")" << std::endl;
197 if (recalculate_peaks_)
200 recalculatePeakBorders_(picked_chroms, best_left, best_right, recalculate_peaks_max_z_);
201 if (peak_apex < best_left || peak_apex > best_right)
204 peak_apex = (best_left + best_right) / 2.0;
207 picked_chroms[chr_idx][peak_idx].setIntensity(0.0);
212 remove_overlapping_features(picked_chroms, best_left, best_right);
215 if (min_peak_width_ > 0.0 && std::fabs(best_right - best_left) < min_peak_width_) {
return mrmFeature;}
217 if (compute_peak_quality_)
220 double qual = computeQuality_(transition_group, picked_chroms, chr_idx, best_left, best_right, outlier);
221 if (qual < min_qual_) {
return mrmFeature; }
230 SpectrumT master_peak_container;
231 const SpectrumT& ref_chromatogram = transition_group.
getChromatograms()[chr_idx];
232 prepareMasterContainer_(ref_chromatogram, master_peak_container, best_left, best_right);
234 double total_intensity = 0;
double total_peak_apices = 0;
double total_xic = 0;
238 for (
typename SpectrumT::const_iterator it = chromatogram.begin(); it != chromatogram.end(); it++)
240 total_xic += it->getIntensity();
244 const SpectrumT used_chromatogram = resampleChromatogram_(chromatogram, master_peak_container, best_left, best_right);
253 double intensity_sum(0.0), rt_sum(0.0);
254 double peak_apex_int = -1;
255 double peak_apex_dist = std::fabs(used_chromatogram.begin()->getMZ() - peak_apex);
257 for (
typename SpectrumT::const_iterator it = used_chromatogram.begin(); it != used_chromatogram.end(); it++)
259 if (it->getMZ() > best_left && it->getMZ() < best_right)
263 p[1] = it->getIntensity();
264 hull_points.push_back(p);
265 if (std::fabs(it->getMZ() - peak_apex) <= peak_apex_dist)
267 peak_apex_int = p[1];
268 peak_apex_dist = std::fabs(it->getMZ() - peak_apex);
270 rt_sum += it->getMZ();
271 intensity_sum += it->getIntensity();
275 if (background_subtraction_ !=
"none")
277 double background = 0;
278 if (background_subtraction_ ==
"smoothed")
294 else if (background_subtraction_ ==
"original")
296 background = calculateBgEstimation_(used_chromatogram, best_left, best_right);
298 intensity_sum -= background;
299 if (intensity_sum < 0)
301 std::cerr <<
"Warning: Intensity was below 0 after background subtraction: " << intensity_sum <<
". Setting it to 0." << std::endl;
306 f.
setRT(picked_chroms[chr_idx][peak_idx].getMZ());
307 f.
setMZ(chromatogram.getMetaValue(
"product_mz"));
312 f.
setMetaValue(
"MZ", chromatogram.getMetaValue(
"product_mz"));
313 f.
setMetaValue(
"native_id", chromatogram.getNativeID());
318 total_intensity += intensity_sum;
319 total_peak_apices += peak_apex_int;
320 mrmFeature.
addFeature(f, chromatogram.getNativeID());
329 const SpectrumT used_chromatogram = resampleChromatogram_(chromatogram, master_peak_container, best_left, best_right);
337 double intensity_sum(0.0), rt_sum(0.0);
338 double peak_apex_int = -1;
339 double peak_apex_dist = std::fabs(used_chromatogram.begin()->getMZ() - peak_apex);
341 for (
typename SpectrumT::const_iterator it = used_chromatogram.begin(); it != used_chromatogram.end(); it++)
343 if (it->getMZ() > best_left && it->getMZ() < best_right)
347 p[1] = it->getIntensity();
348 hull_points.push_back(p);
349 if (std::fabs(it->getMZ() - peak_apex) <= peak_apex_dist)
351 peak_apex_int = p[1];
352 peak_apex_dist = std::fabs(it->getMZ() - peak_apex);
354 rt_sum += it->getMZ();
355 intensity_sum += it->getIntensity();
359 if (chromatogram.metaValueExists(
"precursor_mz"))
361 f.
setMZ(chromatogram.getMetaValue(
"precursor_mz"));
363 f.
setRT(picked_chroms[chr_idx][peak_idx].getMZ());
368 f.
setMetaValue(
"native_id", chromatogram.getNativeID());
374 mrmFeature.
setRT(peak_apex);
380 mrmFeature.
setMetaValue(
"peak_apices_sum", total_peak_apices);
399 template <
typename SpectrumT>
404 for (
Size k = 0;
k < picked_chroms.size();
k++)
406 for (
Size i = 0; i < picked_chroms[
k].size(); i++)
408 if (picked_chroms[
k][i].getMZ() >= best_left && picked_chroms[
k][i].getMZ() <= best_right)
412 picked_chroms[
k][i].setIntensity(0.0);
418 for (
Size k = 0;
k < picked_chroms.size();
k++)
420 for (
Size i = 0; i < picked_chroms[
k].size(); i++)
422 if (picked_chroms[
k][i].getIntensity() <= 0.0) {
continue; }
424 double left = picked_chroms[
k].getFloatDataArrays()[1][i];
425 double right = picked_chroms[
k].getFloatDataArrays()[2][i];
426 if ((left > best_left && left < best_right)
427 || (right > best_left && right < best_right))
431 picked_chroms[
k][i].setIntensity(0.0);
438 void findLargestPeak(std::vector<RichPeakChromatogram>& picked_chroms,
int& chr_idx,
int& peak_idx);
443 void updateMembers_();
464 template <
typename SpectrumT,
typename TransitionT>
466 std::vector<SpectrumT>& picked_chroms,
const int chr_idx,
467 const double best_left,
const double best_right,
String& outlier)
473 double resample_boundary = 15.0;
474 SpectrumT master_peak_container;
475 const SpectrumT& ref_chromatogram = transition_group.
getChromatograms()[chr_idx];
476 prepareMasterContainer_(ref_chromatogram, master_peak_container, best_left - resample_boundary, best_right + resample_boundary);
477 std::vector<std::vector<double> > all_ints;
481 const SpectrumT used_chromatogram = resampleChromatogram_(chromatogram,
482 master_peak_container, best_left - resample_boundary, best_right + resample_boundary);
484 std::vector<double> int_here;
485 for (
Size i = 0; i < used_chromatogram.size(); i++)
487 int_here.push_back(used_chromatogram[i].getIntensity());
489 all_ints.push_back(int_here);
493 std::vector<double> mean_shapes;
494 std::vector<double> mean_coel;
495 for (
Size k = 0;
k < all_ints.size();
k++)
497 std::vector<double> shapes;
498 std::vector<double> coel;
499 for (
Size i = 0; i < all_ints.size(); i++)
501 if (i ==
k) {
continue;}
503 all_ints[
k], all_ints[i], boost::numeric_cast<int>(all_ints[i].size()), 1);
509 shapes.push_back(res_shape);
510 coel.push_back(res_coelution);
516 msc = std::for_each(shapes.begin(), shapes.end(), msc);
517 double shapes_mean = msc.
mean();
518 msc = std::for_each(coel.begin(), coel.end(), msc);
519 double coel_mean = msc.
mean();
523 mean_shapes.push_back(shapes_mean);
524 mean_coel.push_back(coel_mean);
530 int min_index_shape = std::distance(mean_shapes.begin(), std::min_element(mean_shapes.begin(), mean_shapes.end()));
531 int max_index_coel = std::distance(mean_coel.begin(), std::max_element(mean_coel.begin(), mean_coel.end()));
534 int missing_peaks = 0;
535 int multiple_peaks = 0;
538 std::vector<double> left_borders;
539 std::vector<double> right_borders;
540 for (
Size k = 0;
k < picked_chroms.size();
k++)
549 for (
Size i = 0; i < picked_chroms[
k].size(); i++)
551 if (picked_chroms[
k][i].getMZ() >= best_left && picked_chroms[
k][i].getMZ() <= best_right)
554 if (picked_chroms[
k][i].getIntensity() > max_int)
556 max_int = picked_chroms[
k][i].getIntensity() > max_int;
557 l_tmp = picked_chroms[
k].getFloatDataArrays()[1][i];
558 r_tmp = picked_chroms[
k].getFloatDataArrays()[2][i];
563 if (l_tmp > 0.0) left_borders.push_back(l_tmp);
564 if (r_tmp > 0.0) right_borders.push_back(r_tmp);
566 if (pfound == 0) missing_peaks++;
567 if (pfound > 1) multiple_peaks++;
572 LOG_DEBUG <<
" Overall found missing : " << missing_peaks <<
" and multiple : " << multiple_peaks << std::endl;
578 if (min_index_shape == max_index_coel)
580 LOG_DEBUG <<
" element " << min_index_shape <<
" is a candidate for removal ... " << std::endl;
593 double shape_score = std::accumulate(mean_shapes.begin(), mean_shapes.end(), 0.0) / mean_shapes.size();
594 double coel_score = std::accumulate(mean_coel.begin(), mean_coel.end(), 0.0) / mean_coel.size();
595 coel_score = (coel_score - 1.0) / 2.0;
597 double score = shape_score - coel_score - 1.0 * missing_peaks / picked_chroms.size();
599 LOG_DEBUG <<
" computed score " << score <<
" (from " << shape_score <<
600 " - " << coel_score <<
" - " << 1.0 * missing_peaks / picked_chroms.size() <<
")" << std::endl;
614 template <
typename SpectrumT>
621 std::vector<double> left_borders;
622 std::vector<double> right_borders;
623 for (
Size k = 0;
k < picked_chroms.size();
k++)
628 for (
Size i = 0; i < picked_chroms[
k].size(); i++)
630 if (picked_chroms[
k][i].getMZ() >= best_left && picked_chroms[
k][i].getMZ() <= best_right)
632 if (picked_chroms[
k].getFloatDataArrays()[0][i] > max_int)
634 max_int = picked_chroms[
k].getFloatDataArrays()[0][i];
635 left = picked_chroms[
k].getFloatDataArrays()[1][i];
636 right = picked_chroms[
k].getFloatDataArrays()[2][i];
642 left_borders.push_back(left);
643 right_borders.push_back(right);
644 LOG_DEBUG <<
" * " <<
k <<
" left boundary " << left_borders.back() <<
" with int " << max_int << std::endl;
645 LOG_DEBUG <<
" * " <<
k <<
" right boundary " << right_borders.back() <<
" with int " << max_int << std::endl;
650 if (right_borders.empty())
668 mean = std::accumulate(right_borders.begin(), right_borders.end(), 0.0) / (
double) right_borders.size();
669 stdev = std::sqrt(std::inner_product(right_borders.begin(), right_borders.end(), right_borders.begin(), 0.0)
670 / right_borders.size() - mean *
mean);
671 std::sort(right_borders.begin(), right_borders.end());
673 LOG_DEBUG <<
" - Recalculating right peak boundaries " << mean <<
" mean / best "
674 << best_right <<
" std " << stdev <<
" : " << std::fabs(best_right - mean) / stdev
675 <<
" coefficient of variation" << std::endl;
678 if (std::fabs(best_right - mean) / stdev > max_z)
680 best_right = right_borders[right_borders.size() / 2];
681 LOG_DEBUG <<
" - Setting right boundary to " << best_right << std::endl;
685 mean = std::accumulate(left_borders.begin(), left_borders.end(), 0.0) / (
double) left_borders.size();
686 stdev = std::sqrt(std::inner_product(left_borders.begin(), left_borders.end(), left_borders.begin(), 0.0)
687 / left_borders.size() - mean *
mean);
688 std::sort(left_borders.begin(), left_borders.end());
690 LOG_DEBUG <<
" - Recalculating left peak boundaries " << mean <<
" mean / best "
691 << best_left <<
" std " << stdev <<
" : " << std::fabs(best_left - mean) / stdev
692 <<
" coefficient of variation" << std::endl;
695 if (std::fabs(best_left - mean) / stdev > max_z)
697 best_left = left_borders[left_borders.size() / 2];
698 LOG_DEBUG <<
" - Setting left boundary to " << best_left << std::endl;
705 template <
typename SpectrumT>
708 SpectrumT& master_peak_container,
double best_left,
double best_right)
711 typename SpectrumT::const_iterator begin = ref_chromatogram.begin();
712 while (begin != ref_chromatogram.end() && begin->getMZ() < best_left) {begin++; }
713 if (begin != ref_chromatogram.begin()) {begin--; }
715 typename SpectrumT::const_iterator end = begin;
716 while (end != ref_chromatogram.end() && end->getMZ() < best_right) {end++; }
717 if (end != ref_chromatogram.end()) {end++; }
720 master_peak_container.resize(distance(begin, end));
721 typename SpectrumT::iterator it = master_peak_container.begin();
722 for (
typename SpectrumT::const_iterator chrom_it = begin; chrom_it != end; chrom_it++, it++)
724 it->setMZ(chrom_it->getMZ());
729 template <
typename SpectrumT>
731 SpectrumT& master_peak_container,
double best_left,
double best_right)
736 typename SpectrumT::const_iterator begin = chromatogram.begin();
737 while (begin != chromatogram.end() && begin->getMZ() < best_left) {begin++; }
738 if (begin != chromatogram.begin()) {begin--; }
740 typename SpectrumT::const_iterator end = begin;
741 while (end != chromatogram.end() && end->getMZ() < best_right) {end++; }
742 if (end != chromatogram.end()) {end++; }
744 SpectrumT resampled_peak_container = master_peak_container;
746 lresampler.
raster(begin, end, resampled_peak_container.begin(), resampled_peak_container.end());
748 return resampled_peak_container;
759 double calculateBgEstimation_(
const RichPeakChromatogram& chromatogram,
double best_left,
double best_right);
774 #endif // OPENMS_ANALYSIS_OPENSWATH_MRMTRANSITIONGROUPPICKER_H
void raster(SpecT< PeakType > &spectrum)
Applies the resampling algorithm to an MSSpectrum.
Definition: LinearResamplerAlign.h:67
double min_qual_
Definition: MRMTransitionGroupPicker.h:765
String background_subtraction_
Definition: MRMTransitionGroupPicker.h:762
const std::vector< TransitionType > & getTransitions() const
Definition: MRMTransitionGroup.h:125
double recalculate_peaks_max_z_
Definition: MRMTransitionGroupPicker.h:770
A more convenient string class.
Definition: String.h:57
bool compute_peak_quality_
Definition: MRMTransitionGroupPicker.h:764
void setMZ(CoordinateType coordinate)
Mutable access to the m/z coordinate (index 1)
Definition: Peak2D.h:197
void sortByPosition()
Lexicographically sorts the peaks by their position.
Definition: MSSpectrum.h:419
IntensityType getIntensity() const
Definition: Peak2D.h:161
OPENSWATHALGO_DLLAPI XCorrArrayType normalizedCrossCorrelation(std::vector< double > &data1, std::vector< double > &data2, int maxdelay, int lag)
const std::vector< SpectrumType > & getChromatograms() const
Definition: MRMTransitionGroup.h:151
SpectrumType & getPrecursorChromatogram(String key)
Definition: MRMTransitionGroup.h:193
std::vector< PointType > PointArrayType
Definition: ConvexHull2D.h:77
bool isSorted() const
Checks if all peaks are sorted with respect to ascending m/z.
Definition: MSSpectrum.h:481
void prepareMasterContainer_(const SpectrumT &ref_chromatogram, SpectrumT &master_peak_container, double best_left, double best_right)
create an empty master peak container that has the correct mz / RT values set
Definition: MRMTransitionGroupPicker.h:707
The PeakPickerMRM finds peaks a single chromatogram.
Definition: PeakPickerMRM.h:63
int stop_after_feature_
Definition: MRMTransitionGroupPicker.h:767
double mean() const
Definition: StatsHelpers.h:208
Main OpenMS namespace.
Definition: FeatureDeconvolution.h:47
void setParameters(const Param ¶m)
Sets the parameters.
A 2-dimensional hull representation in [counter]clockwise direction - depending on axis labelling...
Definition: ConvexHull2D.h:73
#define LOG_DEBUG
Macro for general debugging information.
Definition: LogStream.h:459
static double mean(IteratorType begin, IteratorType end)
Calculates the mean of a range of values.
Definition: StatisticFunctions.h:131
void setIntensity(IntensityType intensity)
Non-mutable access to the data point intensity (height)
Definition: Peak2D.h:167
void addFeature(Feature &feature, const String &key)
Adds an feature from a single chromatogram into the feature.
double min_peak_width_
Definition: MRMTransitionGroupPicker.h:769
The representation of a 1D spectrum.
Definition: MSSpectrum.h:66
void setRT(CoordinateType coordinate)
Mutable access to the RT coordinate (index 0)
Definition: Peak2D.h:209
The MRMTransitionGroupPicker finds peaks in chromatograms that belong to the same precursors...
Definition: MRMTransitionGroupPicker.h:81
bool recalculate_peaks_
Definition: MRMTransitionGroupPicker.h:763
const std::vector< ConvexHull2D > & getConvexHulls() const
Non-mutable access to the convex hulls.
void sortByIntensity(bool reverse=false)
Lexicographically sorts the peaks by their intensity.
Definition: MSSpectrum.h:342
void setQuality(Size index, QualityType q)
Set the quality in dimension c.
Linear Resampling of raw data with alignment.
Definition: LinearResamplerAlign.h:57
void pickChromatogram(const RichPeakChromatogram &chromatogram, RichPeakChromatogram &picked_chrom)
Finds peaks in a single chromatogram and annotates left/right borders.
void setHullPoints(const PointArrayType &points)
accessor for the outer(!) points (no checking is performed if this is actually a convex hull) ...
double computeQuality_(MRMTransitionGroup< SpectrumT, TransitionT > &transition_group, std::vector< SpectrumT > &picked_chroms, const int chr_idx, const double best_left, const double best_right, String &outlier)
Compute transition group quality (higher score is better)
Definition: MRMTransitionGroupPicker.h:465
An LC-MS feature.
Definition: Feature.h:70
functor to compute the mean and stddev of sequence using the std::foreach algorithm ...
Definition: StatsHelpers.h:170
void addFeature(MRMFeature &feature)
Definition: MRMTransitionGroup.h:213
void addPrecursorFeature(Feature &feature, const String &key)
Adds a precursor feature from a single chromatogram into the feature.
void setOverallQuality(QualityType q)
Set the overall quality.
bool hasPrecursorChromatogram(String key) const
Definition: MRMTransitionGroup.h:198
SpectrumT resampleChromatogram_(const SpectrumT &chromatogram, SpectrumT &master_peak_container, double best_left, double best_right)
use the master container from above to resample a chromatogram at those points stored in the master c...
Definition: MRMTransitionGroupPicker.h:730
void pickTransitionGroup(MRMTransitionGroup< SpectrumT, TransitionT > &transition_group)
Pick a group of chromatograms belonging to the same peptide.
Definition: MRMTransitionGroupPicker.h:115
A base class for all classes handling default parameters.
Definition: DefaultParamHandler.h:92
MRMFeature createMRMFeature(MRMTransitionGroup< SpectrumT, TransitionT > &transition_group, std::vector< SpectrumT > &picked_chroms, int &chr_idx, int &peak_idx)
Create feature from a vector of chromatograms and a specified peak.
Definition: MRMTransitionGroupPicker.h:185
void remove_overlapping_features(std::vector< SpectrumT > &picked_chroms, double best_left, double best_right)
Remove overlapping features.
Definition: MRMTransitionGroupPicker.h:400
OPENSWATHALGO_DLLAPI XCorrArrayType::iterator xcorrArrayGetMaxPeak(XCorrArrayType &array)
Find best peak in an cross-correlation (highest apex)
double stop_after_intensity_ratio_
Definition: MRMTransitionGroupPicker.h:768
A multi-chromatogram MRM feature.
Definition: MRMFeature.h:50
void recalculatePeakBorders_(std::vector< SpectrumT > &picked_chroms, double &best_left, double &best_right, double max_z)
Recalculate the borders of the peak.
Definition: MRMTransitionGroupPicker.h:615
MSSpectrum< ChromatogramPeak > RichPeakChromatogram
Definition: MRMTransitionGroupPicker.h:88
const String & getTransitionGroupID() const
Definition: MRMTransitionGroup.h:115
Not implemented exception.
Definition: Exception.h:437
Size ensureUniqueId()
Assigns a valid unique id, but only if the present one is invalid. Returns 1 if the unique id was cha...
Definition: UniqueIdInterface.h:159