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() )
136 picker_.pickChromatogram(chromatogram, picked_chrom, smoothed_chrom);
138 picked_chroms_.push_back(picked_chrom);
139 smoothed_chroms_.push_back(smoothed_chrom);
147 SpectrumT picked_chrom, smoothed_chrom;
149 String native_id = chromatogram.getNativeID();
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)
break;
180 MRMFeature mrm_feature = createMRMFeature(transition_group, picked_chroms_, smoothed_chroms_, chr_idx, peak_idx);
183 features.push_back(mrm_feature);
187 if (stop_after_feature_ > 0 && cnt > stop_after_feature_) {
break;}
196 for (
Size i = 0; i < features.size(); i++)
200 for (
Size j = 0; j < i; j++)
202 if ((
double)mrm_feature.
getMetaValue(
"leftWidth") >= (
double)features[j].getMetaValue(
"leftWidth") &&
203 (
double)mrm_feature.
getMetaValue(
"rightWidth") <= (
double)features[j].getMetaValue(
"rightWidth") )
215 template <
typename SpectrumT,
typename TransitionT>
217 std::vector<SpectrumT>& picked_chroms,
218 const std::vector<SpectrumT>& smoothed_chroms,
227 double best_left = picked_chroms[chr_idx].getFloatDataArrays()[1][peak_idx];
228 double best_right = picked_chroms[chr_idx].getFloatDataArrays()[2][peak_idx];
229 double peak_apex = picked_chroms[chr_idx][peak_idx].getRT();
230 LOG_DEBUG <<
"**** Creating MRMFeature for peak " << chr_idx <<
" " << peak_idx <<
" " <<
231 picked_chroms[chr_idx][peak_idx] <<
" with borders " << best_left <<
" " <<
232 best_right <<
" (" << best_right - best_left <<
")" << std::endl;
234 if (use_consensus_ && recalculate_peaks_)
237 recalculatePeakBorders_(picked_chroms, best_left, best_right, recalculate_peaks_max_z_);
238 if (peak_apex < best_left || peak_apex > best_right)
241 peak_apex = (best_left + best_right) / 2.0;
245 std::vector< double > left_edges;
246 std::vector< double > right_edges;
247 double min_left = best_left;
248 double max_right = best_right;
254 remove_overlapping_features(picked_chroms, best_left, best_right);
263 for (
Size k = 0;
k < picked_chroms.size();
k++)
265 double peak_apex_dist_min = 1e6;
267 for (
Size i = 0; i < picked_chroms[
k].size(); i++)
270 picked_chroms[
k], picked_chroms[
k].getFloatDataArrays()[1][i], picked_chroms[
k].getFloatDataArrays()[2][i]);
271 if (pa_tmp.
apex_pos > 0.0 && std::fabs(pa_tmp.
apex_pos - peak_apex) < peak_apex_dist_min)
278 double l = best_left;
279 double r = best_right;
282 l = picked_chroms[
k].getFloatDataArrays()[1][min_dist];
283 r = picked_chroms[
k].getFloatDataArrays()[2][min_dist];
284 picked_chroms[
k][min_dist].setIntensity(0.0);
287 left_edges.push_back(l);
288 right_edges.push_back(r);
290 if (l < min_left) {min_left = l;}
291 if (r > max_right) {max_right = r;}
295 picked_chroms[chr_idx][peak_idx].setIntensity(0.0);
298 if (use_consensus_ && min_peak_width_ > 0.0 && std::fabs(best_right - best_left) < min_peak_width_)
303 if (use_consensus_ && compute_peak_quality_)
306 double qual = computeQuality_(transition_group, picked_chroms, chr_idx, best_left, best_right, outlier);
307 if (qual < min_qual_)
320 SpectrumT master_peak_container;
321 const SpectrumT& ref_chromatogram = selectChromHelper_(transition_group, picked_chroms[chr_idx].getNativeID());
322 prepareMasterContainer_(ref_chromatogram, master_peak_container, min_left, max_right);
327 double total_intensity = 0;
double total_peak_apices = 0;
double total_xic = 0;
double total_mi = 0;
331 double local_left = best_left;
332 double local_right = best_right;
335 local_left = left_edges[
k];
336 local_right = right_edges[
k];
339 const SpectrumT& chromatogram = selectChromHelper_(transition_group, transition_group.
getTransitions()[
k].getNativeID());
342 for (
typename SpectrumT::const_iterator it = chromatogram.begin(); it != chromatogram.end(); it++)
344 total_xic += it->getIntensity();
349 double transition_total_xic = 0;
351 for (
typename SpectrumT::const_iterator it = chromatogram.begin(); it != chromatogram.end(); it++)
353 transition_total_xic += it->getIntensity();
357 double transition_total_mi = 0;
358 if (compute_total_mi_)
360 std::vector<double> chrom_vect_id, chrom_vect_det;
361 for (
typename SpectrumT::const_iterator it = chromatogram.begin(); it != chromatogram.end(); it++)
363 chrom_vect_id.push_back(it->getIntensity());
367 int transition_total_mi_norm = 0;
372 const SpectrumT& chromatogram_det = selectChromHelper_(transition_group, transition_group.
getTransitions()[m].getNativeID());
373 chrom_vect_det.clear();
374 for (
typename SpectrumT::const_iterator it = chromatogram_det.begin(); it != chromatogram_det.end(); it++)
376 chrom_vect_det.push_back(it->getIntensity());
379 transition_total_mi_norm++;
382 if (transition_total_mi_norm > 0) { transition_total_mi /= transition_total_mi_norm; }
387 total_mi += transition_total_mi / transition_total_mi_norm;
391 SpectrumT used_chromatogram;
393 if (peak_integration_ ==
"original")
395 used_chromatogram = resampleChromatogram_(chromatogram, master_peak_container, local_left, local_right);
398 else if (peak_integration_ ==
"smoothed" && smoothed_chroms.size() <=
k)
401 "Tried to calculate peak area and height without any smoothed chromatograms");
403 else if (peak_integration_ ==
"smoothed")
405 used_chromatogram = resampleChromatogram_(smoothed_chroms[
k], master_peak_container, local_left, local_right);
410 String(
"Peak integration chromatogram ") + peak_integration_ +
" is not a valid method for MRMTransitionGroupPicker");
419 double peak_integral = pa.
area;
420 double peak_apex_int = pa.
height;
422 if (background_subtraction_ !=
"none")
424 double background{0};
425 double avg_noise_level{0};
426 if ((peak_integration_ ==
"smoothed") && smoothed_chroms.size() <=
k)
429 "Tried to calculate background estimation without any smoothed chromatograms");
431 else if (background_subtraction_ ==
"original")
433 const double intensity_left = chromatogram.PosBegin(local_left)->getIntensity();
434 const double intensity_right = (chromatogram.PosEnd(local_right) - 1)->getIntensity();
435 const UInt n_points = std::distance(chromatogram.PosBegin(local_left), chromatogram.PosEnd(local_right));
436 avg_noise_level = (intensity_right + intensity_left) / 2;
437 background = avg_noise_level * n_points;
439 else if (background_subtraction_ ==
"exact")
442 background = pb.
area;
443 avg_noise_level = pb.
height;
445 peak_integral -= background;
446 peak_apex_int -= avg_noise_level;
447 if (peak_integral < 0) {peak_integral = 0;}
448 if (peak_apex_int < 0) {peak_apex_int = 0;}
451 f.
setMetaValue(
"noise_background_level", avg_noise_level);
454 f.
setRT(picked_chroms[chr_idx][peak_idx].getMZ());
459 if (chromatogram.metaValueExists(
"product_mz"))
461 f.
setMetaValue(
"MZ", chromatogram.getMetaValue(
"product_mz"));
462 f.
setMZ(chromatogram.getMetaValue(
"product_mz"));
466 LOG_WARN <<
"Please set meta value 'product_mz' on chromatogram to populate feature m/z value" << std::endl;
468 f.
setMetaValue(
"native_id", chromatogram.getNativeID());
471 if (compute_total_mi_)
478 total_intensity += peak_integral;
479 total_peak_apices += peak_apex_int;
484 if (compute_peak_shape_metrics_)
505 mrmFeature.
addFeature(f, chromatogram.getNativeID());
515 double local_left = best_left;
516 double local_right = best_right;
517 if (!use_consensus_ && right_edges.size() > prec_idx && left_edges.size() > prec_idx)
519 local_left = left_edges[prec_idx];
520 local_right = right_edges[prec_idx];
523 SpectrumT used_chromatogram;
525 if (peak_integration_ ==
"original")
527 used_chromatogram = resampleChromatogram_(chromatogram, master_peak_container, local_left, local_right);
530 else if (peak_integration_ ==
"smoothed" && smoothed_chroms.size() <= prec_idx)
533 "Tried to calculate peak area and height without any smoothed chromatograms for precursors");
535 else if (peak_integration_ ==
"smoothed")
537 used_chromatogram = resampleChromatogram_(smoothed_chroms[prec_idx], master_peak_container, local_left, local_right);
542 String(
"Peak integration chromatogram ") + peak_integration_ +
" is not a valid method for MRMTransitionGroupPicker");
551 double peak_integral = pa.
area;
552 double peak_apex_int = pa.
height;
554 if (background_subtraction_ !=
"none")
556 double background{0};
557 double avg_noise_level{0};
558 if ((peak_integration_ ==
"smoothed") && smoothed_chroms.size() <= prec_idx)
561 "Tried to calculate background estimation without any smoothed chromatograms");
563 else if (background_subtraction_ ==
"original")
565 const double intensity_left = chromatogram.PosBegin(local_left)->getIntensity();
566 const double intensity_right = (chromatogram.PosEnd(local_right) - 1)->getIntensity();
567 const UInt n_points = std::distance(chromatogram.PosBegin(local_left), chromatogram.PosEnd(local_right));
568 avg_noise_level = (intensity_right + intensity_left) / 2;
569 background = avg_noise_level * n_points;
571 else if (background_subtraction_ ==
"exact")
574 background = pb.
area;
575 avg_noise_level = pb.
height;
577 peak_integral -= background;
578 peak_apex_int -= avg_noise_level;
579 if (peak_integral < 0) {peak_integral = 0;}
580 if (peak_apex_int < 0) {peak_apex_int = 0;}
583 f.
setMetaValue(
"noise_background_level", avg_noise_level);
586 if (chromatogram.metaValueExists(
"precursor_mz"))
588 f.
setMZ(chromatogram.getMetaValue(
"precursor_mz"));
589 mrmFeature.
setMZ(chromatogram.getMetaValue(
"precursor_mz"));
592 f.
setRT(picked_chroms[chr_idx][peak_idx].getMZ());
597 f.
setMetaValue(
"native_id", chromatogram.getNativeID());
602 total_intensity += peak_integral;
608 mrmFeature.
setRT(peak_apex);
614 if (compute_total_mi_)
618 mrmFeature.
setMetaValue(
"peak_apices_sum", total_peak_apices);
637 template <
typename SpectrumT>
642 for (
Size k = 0;
k < picked_chroms.size();
k++)
644 for (
Size i = 0; i < picked_chroms[
k].size(); i++)
646 if (picked_chroms[
k][i].getMZ() >= best_left && picked_chroms[
k][i].getMZ() <= best_right)
650 picked_chroms[
k][i].setIntensity(0.0);
656 for (
Size k = 0;
k < picked_chroms.size();
k++)
658 for (
Size i = 0; i < picked_chroms[
k].size(); i++)
660 if (picked_chroms[
k][i].getIntensity() <= 0.0) {
continue; }
662 double left = picked_chroms[
k].getFloatDataArrays()[1][i];
663 double right = picked_chroms[
k].getFloatDataArrays()[2][i];
664 if ((left > best_left && left < best_right)
665 || (right > best_left && right < best_right))
669 picked_chroms[
k][i].setIntensity(0.0);
676 void findLargestPeak(
const std::vector<MSChromatogram >& picked_chroms,
int& chr_idx,
int& peak_idx);
686 void findWidestPeakIndices(
const std::vector<MSChromatogram>& picked_chroms,
Int& chrom_idx,
Int& point_idx)
const;
691 void updateMembers_()
override;
699 template <
typename SpectrumT,
typename TransitionT>
712 throw Exception::IllegalArgument(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION,
"Did not find chromatogram for id '" + native_id +
"'.");
732 template <
typename SpectrumT,
typename TransitionT>
734 const std::vector<SpectrumT>& picked_chroms,
736 const double best_left,
737 const double best_right,
743 double resample_boundary = resample_boundary_;
744 SpectrumT master_peak_container;
745 const SpectrumT& ref_chromatogram = selectChromHelper_(transition_group, picked_chroms[chr_idx].getNativeID());
746 prepareMasterContainer_(ref_chromatogram, master_peak_container, best_left - resample_boundary, best_right + resample_boundary);
747 std::vector<std::vector<double> > all_ints;
748 for (
Size k = 0;
k < picked_chroms.size();
k++)
750 const SpectrumT chromatogram = selectChromHelper_(transition_group, picked_chroms[
k].getNativeID());
751 const SpectrumT used_chromatogram = resampleChromatogram_(chromatogram,
752 master_peak_container, best_left - resample_boundary, best_right + resample_boundary);
754 std::vector<double> int_here;
755 for (
Size i = 0; i < used_chromatogram.size(); i++)
757 int_here.push_back(used_chromatogram[i].getIntensity());
759 all_ints.push_back(int_here);
763 std::vector<double> mean_shapes;
764 std::vector<double> mean_coel;
765 for (
Size k = 0;
k < all_ints.size();
k++)
767 std::vector<double> shapes;
768 std::vector<double> coel;
769 for (
Size i = 0; i < all_ints.size(); i++)
771 if (i ==
k) {
continue;}
773 all_ints[
k], all_ints[i], boost::numeric_cast<int>(all_ints[i].size()), 1);
779 shapes.push_back(res_shape);
780 coel.push_back(res_coelution);
786 msc = std::for_each(shapes.begin(), shapes.end(), msc);
787 double shapes_mean = msc.
mean();
788 msc = std::for_each(coel.begin(), coel.end(), msc);
789 double coel_mean = msc.
mean();
793 mean_shapes.push_back(shapes_mean);
794 mean_coel.push_back(coel_mean);
800 int min_index_shape = std::distance(mean_shapes.begin(), std::min_element(mean_shapes.begin(), mean_shapes.end()));
801 int max_index_coel = std::distance(mean_coel.begin(), std::max_element(mean_coel.begin(), mean_coel.end()));
804 int missing_peaks = 0;
805 int multiple_peaks = 0;
808 std::vector<double> left_borders;
809 std::vector<double> right_borders;
810 for (
Size k = 0;
k < picked_chroms.size();
k++)
819 for (
Size i = 0; i < picked_chroms[
k].size(); i++)
821 if (picked_chroms[
k][i].getMZ() >= best_left && picked_chroms[
k][i].getMZ() <= best_right)
824 if (picked_chroms[
k][i].getIntensity() > max_int)
826 max_int = picked_chroms[
k][i].getIntensity() > max_int;
827 l_tmp = picked_chroms[
k].getFloatDataArrays()[1][i];
828 r_tmp = picked_chroms[
k].getFloatDataArrays()[2][i];
833 if (l_tmp > 0.0) left_borders.push_back(l_tmp);
834 if (r_tmp > 0.0) right_borders.push_back(r_tmp);
836 if (pfound == 0) missing_peaks++;
837 if (pfound > 1) multiple_peaks++;
842 LOG_DEBUG <<
" Overall found missing : " << missing_peaks <<
" and multiple : " << multiple_peaks << std::endl;
848 if (min_index_shape == max_index_coel)
850 LOG_DEBUG <<
" element " << min_index_shape <<
" is a candidate for removal ... " << std::endl;
851 outlier =
String(picked_chroms[min_index_shape].getNativeID());
863 double shape_score = std::accumulate(mean_shapes.begin(), mean_shapes.end(), 0.0) / mean_shapes.size();
864 double coel_score = std::accumulate(mean_coel.begin(), mean_coel.end(), 0.0) / mean_coel.size();
865 coel_score = (coel_score - 1.0) / 2.0;
867 double score = shape_score - coel_score - 1.0 * missing_peaks / picked_chroms.size();
869 LOG_DEBUG <<
" computed score " << score <<
" (from " << shape_score <<
870 " - " << coel_score <<
" - " << 1.0 * missing_peaks / picked_chroms.size() <<
")" << std::endl;
884 template <
typename SpectrumT>
885 void recalculatePeakBorders_(
const std::vector<SpectrumT>& picked_chroms,
double& best_left,
double& best_right,
double max_z)
891 std::vector<double> left_borders;
892 std::vector<double> right_borders;
893 for (
Size k = 0;
k < picked_chroms.size();
k++)
898 for (
Size i = 0; i < picked_chroms[
k].size(); i++)
900 if (picked_chroms[
k][i].getMZ() >= best_left && picked_chroms[
k][i].getMZ() <= best_right)
902 if (picked_chroms[
k].getFloatDataArrays()[0][i] > max_int)
904 max_int = picked_chroms[
k].getFloatDataArrays()[0][i];
905 left = picked_chroms[
k].getFloatDataArrays()[1][i];
906 right = picked_chroms[
k].getFloatDataArrays()[2][i];
912 left_borders.push_back(left);
913 right_borders.push_back(right);
914 LOG_DEBUG <<
" * " <<
k <<
" left boundary " << left_borders.back() <<
" with int " << max_int << std::endl;
915 LOG_DEBUG <<
" * " <<
k <<
" right boundary " << right_borders.back() <<
" with int " << max_int << std::endl;
920 if (right_borders.empty())
938 mean = std::accumulate(right_borders.begin(), right_borders.end(), 0.0) / (
double) right_borders.size();
939 stdev = std::sqrt(std::inner_product(right_borders.begin(), right_borders.end(), right_borders.begin(), 0.0)
940 / right_borders.size() -
mean *
mean);
941 std::sort(right_borders.begin(), right_borders.end());
943 LOG_DEBUG <<
" - Recalculating right peak boundaries " <<
mean <<
" mean / best "
944 << best_right <<
" std " << stdev <<
" : " << std::fabs(best_right -
mean) / stdev
945 <<
" coefficient of variation" << std::endl;
948 if (std::fabs(best_right -
mean) / stdev > max_z)
950 best_right = right_borders[right_borders.size() / 2];
951 LOG_DEBUG <<
" - Setting right boundary to " << best_right << std::endl;
955 mean = std::accumulate(left_borders.begin(), left_borders.end(), 0.0) / (
double) left_borders.size();
956 stdev = std::sqrt(std::inner_product(left_borders.begin(), left_borders.end(), left_borders.begin(), 0.0)
957 / left_borders.size() -
mean *
mean);
958 std::sort(left_borders.begin(), left_borders.end());
960 LOG_DEBUG <<
" - Recalculating left peak boundaries " <<
mean <<
" mean / best "
961 << best_left <<
" std " << stdev <<
" : " << std::fabs(best_left -
mean) / stdev
962 <<
" coefficient of variation" << std::endl;
965 if (std::fabs(best_left -
mean) / stdev > max_z)
967 best_left = left_borders[left_borders.size() / 2];
968 LOG_DEBUG <<
" - Setting left boundary to " << best_left << std::endl;
990 template <
typename SpectrumT>
992 SpectrumT& master_peak_container,
double left_boundary,
double right_boundary)
999 typename SpectrumT::const_iterator begin = ref_chromatogram.begin();
1000 while (begin != ref_chromatogram.end() && begin->getMZ() < left_boundary) {begin++; }
1001 if (begin != ref_chromatogram.begin()) {begin--; }
1003 typename SpectrumT::const_iterator end = begin;
1004 while (end != ref_chromatogram.end() && end->getMZ() < right_boundary) {end++; }
1005 if (end != ref_chromatogram.end()) {end++; }
1008 master_peak_container.resize(distance(begin, end));
1009 typename SpectrumT::iterator it = master_peak_container.begin();
1010 for (
typename SpectrumT::const_iterator chrom_it = begin; chrom_it != end; chrom_it++, it++)
1012 it->setMZ(chrom_it->getMZ());
1026 template <
typename SpectrumT>
1028 const SpectrumT& master_peak_container,
double left_boundary,
double right_boundary)
1033 typename SpectrumT::const_iterator begin = chromatogram.begin();
1034 while (begin != chromatogram.end() && begin->getMZ() < left_boundary) {begin++;}
1035 if (begin != chromatogram.begin()) {begin--;}
1037 typename SpectrumT::const_iterator end = begin;
1038 while (end != chromatogram.end() && end->getMZ() < right_boundary) {end++;}
1039 if (end != chromatogram.end()) {end++;}
1041 SpectrumT resampled_peak_container = master_peak_container;
1043 lresampler.
raster(begin, end, resampled_peak_container.begin(), resampled_peak_container.end());
1045 return resampled_peak_container;