112 template <
typename SpectrumT,
typename TransitionT>
118 std::vector<MSChromatogram > picked_chroms;
119 std::vector<MSChromatogram > smoothed_chroms;
130 !transition_group.
getTransition(native_id).isDetectingTransition() )
137 picker_.pickChromatogram(chromatogram, picked_chrom, smoothed_chrom);
139 picked_chroms.push_back(std::move(picked_chrom));
140 smoothed_chroms.push_back(std::move(smoothed_chrom));
148 SpectrumT picked_chrom, smoothed_chrom;
151 picker_.pickChromatogram(chromatogram, picked_chrom, smoothed_chrom);
152 picked_chrom.sortByIntensity();
153 picked_chroms.push_back(picked_chrom);
154 smoothed_chroms.push_back(smoothed_chrom);
162 int chr_idx, peak_idx, cnt = 0;
163 std::vector<MRMFeature> features;
166 chr_idx = -1; peak_idx = -1;
168 if (boundary_selection_method_ ==
"largest")
170 findLargestPeak(picked_chroms, chr_idx, peak_idx);
172 else if (boundary_selection_method_ ==
"widest")
174 findWidestPeakIndices(picked_chroms, chr_idx, peak_idx);
177 if (chr_idx == -1 && peak_idx == -1)
179 OPENMS_LOG_DEBUG <<
"**** MRMTransitionGroupPicker : no more peaks left" << picked_chroms.size() << std::endl;
184 MRMFeature mrm_feature = createMRMFeature(transition_group, picked_chroms, smoothed_chroms, chr_idx, peak_idx);
185 double total_xic = 0;
190 features.push_back(std::move(mrm_feature));
194 if (stop_after_feature_ > 0 && cnt > stop_after_feature_) {
break;}
195 if (intensity > 0 && intensity / total_xic < stop_after_intensity_ratio_)
202 for (
Size i = 0; i < features.size(); i++)
206 for (
Size j = 0; j < i; j++)
208 if ((
double)mrm_feature.
getMetaValue(
"leftWidth") >= (
double)features[j].getMetaValue(
"leftWidth") &&
221 template <
typename SpectrumT,
typename TransitionT>
223 std::vector<SpectrumT>& picked_chroms,
224 const std::vector<SpectrumT>& smoothed_chroms,
235 double peak_apex = picked_chroms[chr_idx][peak_idx].getRT();
236 OPENMS_LOG_DEBUG <<
"**** Creating MRMFeature for peak " << chr_idx <<
" " << peak_idx <<
" " <<
237 picked_chroms[chr_idx][peak_idx] <<
" with borders " << best_left <<
" " <<
238 best_right <<
" (" << best_right - best_left <<
")" << std::endl;
240 if (use_consensus_ && recalculate_peaks_)
243 recalculatePeakBorders_(picked_chroms, best_left, best_right, recalculate_peaks_max_z_);
244 if (peak_apex < best_left || peak_apex > best_right)
247 peak_apex = (best_left + best_right) / 2.0;
251 std::vector< double > left_edges;
252 std::vector< double > right_edges;
253 double min_left = best_left;
254 double max_right = best_right;
260 remove_overlapping_features(picked_chroms, best_left, best_right);
264 pickApex(picked_chroms, best_left, best_right, peak_apex,
265 min_left, max_right, left_edges, right_edges);
268 picked_chroms[chr_idx][peak_idx].setIntensity(0.0);
273 if (min_peak_width_ > 0.0 && std::fabs(best_right - best_left) < min_peak_width_)
278 if (compute_peak_quality_)
281 double qual = computeQuality_(transition_group, picked_chroms, chr_idx, best_left, best_right, outlier);
282 if (qual < min_qual_)
286 mrmFeature.setMetaValue(
"potentialOutlier", outlier);
287 mrmFeature.setMetaValue(
"initialPeakQuality", qual);
288 mrmFeature.setOverallQuality(qual);
296 SpectrumT master_peak_container;
297 const SpectrumT& ref_chromatogram = selectChromHelper_(transition_group, picked_chroms[chr_idx].getNativeID());
298 prepareMasterContainer_(ref_chromatogram, master_peak_container, min_left, max_right);
303 double total_intensity = 0;
double total_peak_apices = 0;
double total_xic = 0;
double total_mi = 0;
304 pickFragmentChromatograms(transition_group, picked_chroms, mrmFeature, smoothed_chroms,
305 best_left, best_right, use_consensus_,
306 total_intensity, total_xic, total_mi, total_peak_apices,
307 master_peak_container, left_edges, right_edges,
312 pickPrecursorChromatograms(transition_group,
313 picked_chroms, mrmFeature, smoothed_chroms,
314 best_left, best_right, use_consensus_,
315 total_intensity, master_peak_container, left_edges, right_edges,
318 mrmFeature.setRT(peak_apex);
319 mrmFeature.setIntensity(total_intensity);
321 mrmFeature.setMetaValue(
"leftWidth", best_left);
322 mrmFeature.setMetaValue(
"rightWidth", best_right);
323 mrmFeature.setMetaValue(
"total_xic", total_xic);
324 if (compute_total_mi_)
326 mrmFeature.setMetaValue(
"total_mi", total_mi);
328 mrmFeature.setMetaValue(
"peak_apices_sum", total_peak_apices);
330 mrmFeature.ensureUniqueId();
346 template <
typename SpectrumT>
347 void pickApex(std::vector<SpectrumT>& picked_chroms,
348 const double best_left,
const double best_right,
const double peak_apex,
349 double &min_left,
double &max_right,
350 std::vector< double > & left_edges, std::vector< double > & right_edges)
352 for (
Size k = 0;
k < picked_chroms.size();
k++)
354 double peak_apex_dist_min = std::numeric_limits<double>::max();
356 for (
Size i = 0; i < picked_chroms[
k].size(); i++)
362 if (pa_tmp.
apex_pos > 0.0 && std::fabs(pa_tmp.
apex_pos - peak_apex) < peak_apex_dist_min)
364 peak_apex_dist_min = std::fabs(pa_tmp.
apex_pos - peak_apex);
370 double l = best_left;
371 double r = best_right;
376 picked_chroms[
k][min_dist].setIntensity(0.0);
379 left_edges.push_back(l);
380 right_edges.push_back(r);
382 if (l < min_left) {min_left = l;}
383 if (r > max_right) {max_right = r;}
387 template <
typename SpectrumT,
typename TransitionT>
389 const std::vector<SpectrumT>& picked_chroms,
391 const std::vector<SpectrumT>& smoothed_chroms,
392 const double best_left,
const double best_right,
393 const bool use_consensus_,
394 double & total_intensity,
397 double & total_peak_apices,
398 const SpectrumT & master_peak_container,
399 const std::vector< double > & left_edges,
400 const std::vector< double > & right_edges,
407 double local_left = best_left;
408 double local_right = best_right;
417 "When using non-censensus peak picker, all transitions need to be detecting transitions.");
419 local_left = left_edges[
k];
420 local_right = right_edges[
k];
423 const SpectrumT& chromatogram = selectChromHelper_(transition_group, transition_group.
getTransitions()[
k].getNativeID());
426 for (
typename SpectrumT::const_iterator it = chromatogram.begin(); it != chromatogram.end(); it++)
428 total_xic += it->getIntensity();
433 double transition_total_xic = 0;
435 for (
typename SpectrumT::const_iterator it = chromatogram.begin(); it != chromatogram.end(); it++)
437 transition_total_xic += it->getIntensity();
441 double transition_total_mi = 0;
442 if (compute_total_mi_)
444 std::vector<double> chrom_vect_id, chrom_vect_det;
445 for (
typename SpectrumT::const_iterator it = chromatogram.begin(); it != chromatogram.end(); it++)
447 chrom_vect_id.push_back(it->getIntensity());
451 int transition_total_mi_norm = 0;
456 const SpectrumT& chromatogram_det = selectChromHelper_(transition_group, transition_group.
getTransitions()[m].getNativeID());
457 chrom_vect_det.clear();
458 for (
typename SpectrumT::const_iterator it = chromatogram_det.begin(); it != chromatogram_det.end(); it++)
460 chrom_vect_det.push_back(it->getIntensity());
463 transition_total_mi_norm++;
466 if (transition_total_mi_norm > 0) { transition_total_mi /= transition_total_mi_norm; }
471 total_mi += transition_total_mi / transition_total_mi_norm;
475 SpectrumT used_chromatogram;
477 if (peak_integration_ ==
"original")
479 used_chromatogram = resampleChromatogram_(chromatogram, master_peak_container, local_left, local_right);
481 else if (peak_integration_ ==
"smoothed")
483 if (smoothed_chroms.size() <=
k)
486 "Tried to calculate peak area and height without any smoothed chromatograms");
488 used_chromatogram = resampleChromatogram_(smoothed_chroms[
k], master_peak_container, local_left, local_right);
493 String(
"Peak integration chromatogram ") + peak_integration_ +
" is not a valid method for MRMTransitionGroupPicker");
502 double peak_integral = pa.
area;
503 double peak_apex_int = pa.
height;
505 if (background_subtraction_ !=
"none")
507 double background{0};
508 double avg_noise_level{0};
509 if (background_subtraction_ ==
"original")
511 const double intensity_left = chromatogram.PosBegin(local_left)->getIntensity();
512 const double intensity_right = (chromatogram.PosEnd(local_right) - 1)->getIntensity();
513 const UInt n_points = std::distance(chromatogram.PosBegin(local_left), chromatogram.PosEnd(local_right));
514 avg_noise_level = (intensity_right + intensity_left) / 2;
515 background = avg_noise_level * n_points;
517 else if (background_subtraction_ ==
"exact")
520 background = pb.
area;
521 avg_noise_level = pb.
height;
523 peak_integral -= background;
524 peak_apex_int -= avg_noise_level;
525 if (peak_integral < 0) {peak_integral = 0;}
526 if (peak_apex_int < 0) {peak_apex_int = 0;}
529 f.
setMetaValue(
"noise_background_level", avg_noise_level);
532 f.
setRT(picked_chroms[chr_idx][peak_idx].getMZ());
538 f.
setMZ(chromatogram.getProduct().getMZ());
539 mrmFeature.
setMZ(chromatogram.getPrecursor().getMZ());
541 if (chromatogram.metaValueExists(
"product_mz"))
543 f.
setMetaValue(
"MZ", chromatogram.getMetaValue(
"product_mz"));
544 f.
setMZ(chromatogram.getMetaValue(
"product_mz"));
547 f.
setMetaValue(
"native_id", chromatogram.getNativeID());
550 if (compute_total_mi_)
557 total_intensity += peak_integral;
558 total_peak_apices += peak_apex_int;
565 if (compute_peak_shape_metrics_)
584 mrmFeature.
addFeature(f, chromatogram.getNativeID());
588 template <
typename SpectrumT,
typename TransitionT>
590 const std::vector<SpectrumT>& picked_chroms,
592 const std::vector<SpectrumT>& smoothed_chroms,
593 const double best_left,
const double best_right,
594 const bool use_consensus_,
595 double & total_intensity,
596 const SpectrumT & master_peak_container,
597 const std::vector< double > & left_edges,
598 const std::vector< double > & right_edges,
610 double local_left = best_left;
611 double local_right = best_right;
612 if (!use_consensus_ && right_edges.size() > prec_idx && left_edges.size() > prec_idx)
614 local_left = left_edges[prec_idx];
615 local_right = right_edges[prec_idx];
618 SpectrumT used_chromatogram;
620 if (peak_integration_ ==
"original")
622 used_chromatogram = resampleChromatogram_(chromatogram, master_peak_container, local_left, local_right);
625 else if (peak_integration_ ==
"smoothed" && smoothed_chroms.size() <= prec_idx)
628 "Tried to calculate peak area and height without any smoothed chromatograms for precursors");
630 else if (peak_integration_ ==
"smoothed")
632 used_chromatogram = resampleChromatogram_(smoothed_chroms[prec_idx], master_peak_container, local_left, local_right);
637 String(
"Peak integration chromatogram ") + peak_integration_ +
" is not a valid method for MRMTransitionGroupPicker");
646 double peak_integral = pa.
area;
647 double peak_apex_int = pa.
height;
649 if (background_subtraction_ !=
"none")
651 double background{0};
652 double avg_noise_level{0};
653 if ((peak_integration_ ==
"smoothed") && smoothed_chroms.size() <= prec_idx)
656 "Tried to calculate background estimation without any smoothed chromatograms");
658 else if (background_subtraction_ ==
"original")
660 const double intensity_left = chromatogram.PosBegin(local_left)->getIntensity();
661 const double intensity_right = (chromatogram.PosEnd(local_right) - 1)->getIntensity();
662 const UInt n_points = std::distance(chromatogram.PosBegin(local_left), chromatogram.PosEnd(local_right));
663 avg_noise_level = (intensity_right + intensity_left) / 2;
664 background = avg_noise_level * n_points;
666 else if (background_subtraction_ ==
"exact")
669 background = pb.
area;
670 avg_noise_level = pb.
height;
672 peak_integral -= background;
673 peak_apex_int -= avg_noise_level;
674 if (peak_integral < 0) {peak_integral = 0;}
675 if (peak_apex_int < 0) {peak_apex_int = 0;}
678 f.
setMetaValue(
"noise_background_level", avg_noise_level);
681 f.
setMZ(chromatogram.getPrecursor().getMZ());
682 if (
k == 0) {mrmFeature.
setMZ(chromatogram.getPrecursor().getMZ());}
684 if (chromatogram.metaValueExists(
"precursor_mz"))
686 f.
setMZ(chromatogram.getMetaValue(
"precursor_mz"));
687 if (
k == 0) {mrmFeature.
setMZ(chromatogram.getMetaValue(
"precursor_mz"));}
690 f.
setRT(picked_chroms[chr_idx][peak_idx].getMZ());
695 f.
setMetaValue(
"native_id", chromatogram.getNativeID());
700 total_intensity += peak_integral;
720 template <
typename SpectrumT>
724 for (
Size k = 0;
k < picked_chroms.size();
k++)
726 for (
Size i = 0; i < picked_chroms[
k].size(); i++)
728 if (picked_chroms[
k][i].getMZ() >= best_left && picked_chroms[
k][i].getMZ() <= best_right)
730 picked_chroms[
k][i].setIntensity(0.0);
736 for (
Size k = 0;
k < picked_chroms.size();
k++)
738 for (
Size i = 0; i < picked_chroms[
k].size(); i++)
740 if (picked_chroms[
k][i].getIntensity() <= 0.0) {
continue; }
744 if ((left > best_left && left < best_right)
745 || (right > best_left && right < best_right))
747 picked_chroms[
k][i].setIntensity(0.0);
754 void findLargestPeak(
const std::vector<MSChromatogram >& picked_chroms,
int& chr_idx,
int& peak_idx);
764 void findWidestPeakIndices(
const std::vector<MSChromatogram>& picked_chroms,
Int& chrom_idx,
Int& point_idx)
const;
769 void updateMembers_()
override;
777 template <
typename SpectrumT,
typename TransitionT>
790 throw Exception::IllegalArgument(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION,
"Did not find chromatogram for id '" + native_id +
"'.");
810 template <
typename SpectrumT,
typename TransitionT>
812 const std::vector<SpectrumT>& picked_chroms,
814 const double best_left,
815 const double best_right,
821 double resample_boundary = resample_boundary_;
822 SpectrumT master_peak_container;
823 const SpectrumT& ref_chromatogram = selectChromHelper_(transition_group, picked_chroms[chr_idx].getNativeID());
824 prepareMasterContainer_(ref_chromatogram, master_peak_container, best_left - resample_boundary, best_right + resample_boundary);
825 std::vector<std::vector<double> > all_ints;
826 for (
Size k = 0;
k < picked_chroms.size();
k++)
828 const SpectrumT& chromatogram = selectChromHelper_(transition_group, picked_chroms[
k].getNativeID());
829 const SpectrumT used_chromatogram = resampleChromatogram_(chromatogram,
830 master_peak_container, best_left - resample_boundary, best_right + resample_boundary);
832 std::vector<double> int_here;
833 for (
const auto& peak : used_chromatogram) int_here.push_back(peak.getIntensity());
835 double tic = std::accumulate(int_here.begin(), int_here.end(), 0.0);
836 if (tic > 0.0) all_ints.push_back(int_here);
840 std::vector<double> mean_shapes;
841 std::vector<double> mean_coel;
842 for (
Size k = 0;
k < all_ints.size();
k++)
844 std::vector<double> shapes;
845 std::vector<double> coel;
846 for (
Size i = 0; i < all_ints.size(); i++)
848 if (i ==
k) {
continue;}
850 all_ints[
k], all_ints[i], boost::numeric_cast<int>(all_ints[i].size()), 1);
856 shapes.push_back(res_shape);
857 coel.push_back(res_coelution);
863 msc = std::for_each(shapes.begin(), shapes.end(), msc);
864 double shapes_mean = msc.
mean();
865 msc = std::for_each(coel.begin(), coel.end(), msc);
866 double coel_mean = msc.
mean();
870 mean_shapes.push_back(shapes_mean);
871 mean_coel.push_back(coel_mean);
877 int min_index_shape = std::distance(mean_shapes.begin(), std::min_element(mean_shapes.begin(), mean_shapes.end()));
878 int max_index_coel = std::distance(mean_coel.begin(), std::max_element(mean_coel.begin(), mean_coel.end()));
881 int missing_peaks = 0;
882 int multiple_peaks = 0;
885 std::vector<double> left_borders;
886 std::vector<double> right_borders;
887 for (
Size k = 0;
k < picked_chroms.size();
k++)
896 for (
Size i = 0; i < picked_chroms[
k].size(); i++)
898 if (picked_chroms[
k][i].getMZ() >= best_left && picked_chroms[
k][i].getMZ() <= best_right)
901 if (picked_chroms[
k][i].getIntensity() > max_int)
903 max_int = picked_chroms[
k][i].getIntensity() > max_int;
910 if (l_tmp > 0.0) left_borders.push_back(l_tmp);
911 if (r_tmp > 0.0) right_borders.push_back(r_tmp);
913 if (pfound == 0) missing_peaks++;
914 if (pfound > 1) multiple_peaks++;
919 OPENMS_LOG_DEBUG <<
" Overall found missing : " << missing_peaks <<
" and multiple : " << multiple_peaks << std::endl;
925 if (min_index_shape == max_index_coel)
927 OPENMS_LOG_DEBUG <<
" element " << min_index_shape <<
" is a candidate for removal ... " << std::endl;
928 outlier =
String(picked_chroms[min_index_shape].getNativeID());
940 double shape_score = std::accumulate(mean_shapes.begin(), mean_shapes.end(), 0.0) / mean_shapes.size();
941 double coel_score = std::accumulate(mean_coel.begin(), mean_coel.end(), 0.0) / mean_coel.size();
942 coel_score = (coel_score - 1.0) / 2.0;
944 double score = shape_score - coel_score - 1.0 * missing_peaks / picked_chroms.size();
946 OPENMS_LOG_DEBUG <<
" computed score " << score <<
" (from " << shape_score <<
947 " - " << coel_score <<
" - " << 1.0 * missing_peaks / picked_chroms.size() <<
")" << std::endl;
961 template <
typename SpectrumT>
962 void recalculatePeakBorders_(
const std::vector<SpectrumT>& picked_chroms,
double& best_left,
double& best_right,
double max_z)
968 std::vector<double> left_borders;
969 std::vector<double> right_borders;
970 for (
Size k = 0;
k < picked_chroms.size();
k++)
975 for (
Size i = 0; i < picked_chroms[
k].size(); i++)
977 if (picked_chroms[
k][i].getMZ() >= best_left && picked_chroms[
k][i].getMZ() <= best_right)
989 left_borders.push_back(left);
990 right_borders.push_back(right);
991 OPENMS_LOG_DEBUG <<
" * " <<
k <<
" left boundary " << left_borders.back() <<
" with int " << max_int << std::endl;
992 OPENMS_LOG_DEBUG <<
" * " <<
k <<
" right boundary " << right_borders.back() <<
" with int " << max_int << std::endl;
997 if (right_borders.empty())
1015 mean = std::accumulate(right_borders.begin(), right_borders.end(), 0.0) / (
double) right_borders.size();
1016 stdev = std::sqrt(std::inner_product(right_borders.begin(), right_borders.end(), right_borders.begin(), 0.0)
1017 / right_borders.size() -
mean *
mean);
1018 std::sort(right_borders.begin(), right_borders.end());
1021 << best_right <<
" std " << stdev <<
" : " << std::fabs(best_right -
mean) / stdev
1022 <<
" coefficient of variation" << std::endl;
1025 if (std::fabs(best_right -
mean) / stdev > max_z)
1027 best_right = right_borders[right_borders.size() / 2];
1028 OPENMS_LOG_DEBUG <<
" - Setting right boundary to " << best_right << std::endl;
1032 mean = std::accumulate(left_borders.begin(), left_borders.end(), 0.0) / (
double) left_borders.size();
1033 stdev = std::sqrt(std::inner_product(left_borders.begin(), left_borders.end(), left_borders.begin(), 0.0)
1034 / left_borders.size() -
mean *
mean);
1035 std::sort(left_borders.begin(), left_borders.end());
1038 << best_left <<
" std " << stdev <<
" : " << std::fabs(best_left -
mean) / stdev
1039 <<
" coefficient of variation" << std::endl;
1042 if (std::fabs(best_left -
mean) / stdev > max_z)
1044 best_left = left_borders[left_borders.size() / 2];
1045 OPENMS_LOG_DEBUG <<
" - Setting left boundary to " << best_left << std::endl;
1067 template <
typename SpectrumT>
1069 SpectrumT& master_peak_container,
double left_boundary,
double right_boundary)
1076 typename SpectrumT::const_iterator begin = ref_chromatogram.begin();
1077 while (begin != ref_chromatogram.end() && begin->getMZ() < left_boundary) {begin++; }
1078 if (begin != ref_chromatogram.begin()) {begin--; }
1080 typename SpectrumT::const_iterator end = begin;
1081 while (end != ref_chromatogram.end() && end->getMZ() < right_boundary) {end++; }
1082 if (end != ref_chromatogram.end()) {end++; }
1085 master_peak_container.resize(distance(begin, end));
1086 typename SpectrumT::iterator it = master_peak_container.begin();
1087 for (
typename SpectrumT::const_iterator chrom_it = begin; chrom_it != end; chrom_it++, it++)
1089 it->setMZ(chrom_it->getMZ());
1103 template <
typename SpectrumT>
1105 const SpectrumT& master_peak_container,
double left_boundary,
double right_boundary)
1110 typename SpectrumT::const_iterator begin = chromatogram.begin();
1111 while (begin != chromatogram.end() && begin->getMZ() < left_boundary) {begin++;}
1112 if (begin != chromatogram.begin()) {begin--;}
1114 typename SpectrumT::const_iterator end = begin;
1115 while (end != chromatogram.end() && end->getMZ() < right_boundary) {end++;}
1116 if (end != chromatogram.end()) {end++;}
1118 SpectrumT resampled_peak_container = master_peak_container;
1120 lresampler.
raster(begin, end, resampled_peak_container.begin(), resampled_peak_container.end());
1122 return resampled_peak_container;