83 defaults_.setValue(
"rt_tolerance", 10.0,
"Maximal RT distance (in [s]) for two spectra's precursors.");
84 defaults_.setValue(
"mz_tolerance", 1.0,
"Maximal m/z distance (in Da) for two spectra's precursors.");
90 rt_max_ = (
double) param_.getValue(
"rt_tolerance");
91 mz_max_ = (
double) param_.getValue(
"mz_tolerance");
99 return 1 - ((d_rt / rt_max_ + d_mz / mz_max_) / 2);
106 double d_rt = fabs(first.
getRT() - second.
getRT());
107 double d_mz = fabs(first.
getMZ() - second.
getMZ());
109 if (d_rt > rt_max_ || d_mz > mz_max_)
115 double sim = getSimilarity(d_rt, d_mz);
155 template <
typename MapType>
158 IntList ms_levels = param_.getValue(
"block_method:ms_levels");
159 Int rt_block_size(param_.getValue(
"block_method:rt_block_size"));
160 double rt_max_length = (param_.getValue(
"block_method:rt_max_length"));
162 if (rt_max_length == 0)
164 rt_max_length = (std::numeric_limits<double>::max)();
167 for (IntList::iterator it_mslevel = ms_levels.begin(); it_mslevel < ms_levels.end(); ++it_mslevel)
171 SignedSize block_size_count(rt_block_size + 1);
172 Size idx_spectrum(0);
175 if (
Int(it1->getMSLevel()) == *it_mslevel)
178 if (++block_size_count >= rt_block_size ||
179 exp[idx_spectrum].getRT() - exp[idx_block].getRT() > rt_max_length)
181 block_size_count = 0;
182 idx_block = idx_spectrum;
186 spectra_to_merge[idx_block].push_back(idx_spectrum);
193 if (block_size_count == 0)
195 spectra_to_merge[idx_block] = std::vector<Size>();
199 mergeSpectra_(exp, spectra_to_merge, *it_mslevel);
208 template <
typename MapType>
214 std::vector<BinaryTreeNode> tree;
218 std::vector<BaseFeature> data;
220 for (
Size i = 0; i < exp.
size(); ++i)
222 if (exp[i].getMSLevel() != 2)
228 index_mapping[data.size()] = i;
232 bf.
setRT(exp[i].getRT());
233 std::vector<Precursor> pcs = exp[i].getPrecursors();
240 LOG_WARN <<
"More than one precursor found. Using first one!" << std::endl;
242 bf.
setMZ(pcs[0].getMZ());
245 data_size = data.size();
260 std::vector<std::vector<Size> > clusters;
263 for (
Size ii = 0; ii < tree.size(); ++ii)
265 if (tree[ii].distance >= 1)
267 tree[ii].distance = -1;
269 if (tree[ii].distance != -1)
274 ca.
cut(data_size - node_count, tree, clusters);
282 for (
Size i_outer = 0; i_outer < clusters.size(); ++i_outer)
284 if (clusters[i_outer].size() <= 1)
289 Size cl_index0 = clusters[i_outer][0];
290 spectra_to_merge[index_mapping[cl_index0]] = std::vector<Size>();
292 for (
Size i_inner = 1; i_inner < clusters[i_outer].size(); ++i_inner)
294 Size cl_index = clusters[i_outer][i_inner];
295 spectra_to_merge[index_mapping[cl_index0]].push_back(index_mapping[cl_index]);
300 mergeSpectra_(exp, spectra_to_merge, 2);
313 template <
typename MapType>
317 int ms_level = param_.getValue(
"average_gaussian:ms_level");
318 if (average_type ==
"tophat")
320 ms_level = param_.getValue(
"average_tophat:ms_level");
324 String spectrum_type = param_.getValue(
"average_gaussian:spectrum_type");
325 if (average_type ==
"tophat")
327 spectrum_type = param_.getValue(
"average_tophat:spectrum_type");
331 double fwhm(param_.getValue(
"average_gaussian:rt_FWHM"));
332 double factor = -4 * log(2.0) / (fwhm * fwhm);
333 double cutoff(param_.getValue(
"average_gaussian:cutoff"));
336 bool unit(param_.getValue(
"average_tophat:rt_unit") ==
"scans");
337 double range(param_.getValue(
"average_tophat:rt_range"));
338 double range_seconds = range / 2;
339 int range_scans = range;
340 if ((range_scans % 2) == 0)
344 range_scans = (range_scans - 1) / 2;
352 if (
Int(it_rt->getMSLevel()) == ms_level)
363 terminate_now =
false;
364 while (it_rt_2 != exp.
end() && !terminate_now)
366 if (
Int(it_rt_2->getMSLevel()) == ms_level)
369 if (average_type ==
"gaussian")
371 weight = std::exp(factor * pow(it_rt_2->getRT() - it_rt->getRT(), 2));
373 std::pair<Size, double> p(m, weight);
374 spectra_to_average_over[n].push_back(p);
377 if (average_type ==
"gaussian")
380 terminate_now = std::exp(factor * pow(it_rt_2->getRT() - it_rt->getRT(), 2)) < cutoff;
385 terminate_now = (steps > range_scans);
390 terminate_now = (std::abs(it_rt_2->getRT() - it_rt->getRT()) > range_seconds);
400 terminate_now =
false;
401 while (it_rt_2 != exp.
begin() && !terminate_now)
403 if (
Int(it_rt_2->getMSLevel()) == ms_level)
406 if (average_type ==
"gaussian")
408 weight = std::exp(factor * pow(it_rt_2->getRT() - it_rt->getRT(), 2));
410 std::pair<Size, double> p(m, weight);
411 spectra_to_average_over[n].push_back(p);
414 if (average_type ==
"gaussian")
417 terminate_now = std::exp(factor * pow(it_rt_2->getRT() - it_rt->getRT(), 2)) < cutoff;
422 terminate_now = (steps > range_scans);
427 terminate_now = (std::abs(it_rt_2->getRT() - it_rt->getRT()) > range_seconds);
441 for (std::vector<std::pair<Size, double> >::const_iterator it2 = it->second.begin(); it2 != it->second.end(); ++it2)
446 for (std::vector<std::pair<Size, double> >::iterator it2 = it->second.begin(); it2 != it->second.end(); ++it2)
448 (*it2).second /=
sum;
454 if (spectrum_type ==
"automatic")
456 Size idx = spectra_to_average_over.begin()->first;
457 type = exp[idx].getType();
463 else if (spectrum_type ==
"profile")
467 else if (spectrum_type ==
"centroid")
475 averageCentroidSpectra_(exp, spectra_to_average_over, ms_level);
479 averageProfileSpectra_(exp, spectra_to_average_over, ms_level);
501 template <
typename MapType>
504 double mz_binning_width(param_.getValue(
"mz_binning_width"));
505 String mz_binning_unit(param_.getValue(
"mz_binning_width_unit"));
511 std::set<Size> merged_indices;
516 p.
setValue(
"tolerance", mz_binning_width);
517 if (!(mz_binning_unit ==
"Da" || mz_binning_unit ==
"ppm"))
522 p.
setValue(
"is_relative_tolerance", mz_binning_unit ==
"Da" ?
"false" :
"true");
524 std::vector<std::pair<Size, Size> > alignment;
526 Size count_peaks_aligned(0);
527 Size count_peaks_overall(0);
530 for (
auto it = spectra_to_merge.begin(); it != spectra_to_merge.end(); ++it)
532 ++cluster_sizes[it->second.size() + 1];
538 merged_indices.insert(it->first);
541 double rt_average = consensus_spec.
getRT();
542 double precursor_mz_average = 0.0;
543 Size precursor_count(0);
546 precursor_mz_average = consensus_spec.
getPrecursors()[0].getMZ();
550 count_peaks_overall += consensus_spec.size();
553 for (
auto sit = it->second.begin(); sit != it->second.end(); ++sit)
555 consensus_spec.
unify(exp[*sit]);
556 merged_indices.insert(*sit);
558 rt_average += exp[*sit].getRT();
559 if (ms_level >= 2 && exp[*sit].getPrecursors().size() > 0)
561 precursor_mz_average += exp[*sit].getPrecursors()[0].getMZ();
568 count_peaks_aligned += alignment.size();
569 count_peaks_overall += exp[*sit].
size();
572 Size spec_b_index(0);
575 Size spec_a = consensus_spec.size(), spec_b = exp[*sit].
size(), align_size = alignment.size();
576 for (
auto pit = exp[*sit].begin(); pit != exp[*sit].
end(); ++pit)
578 if (alignment.size() == 0 || alignment[align_index].second != spec_b_index)
581 consensus_spec.push_back(*pit);
587 Size copy_of_align_index(align_index);
589 while (alignment.size() > 0 &&
590 copy_of_align_index < alignment.size() &&
591 alignment[copy_of_align_index].second == spec_b_index)
593 ++copy_of_align_index;
597 while (alignment.size() > 0 &&
598 align_index < alignment.size() &&
599 alignment[align_index].second == spec_b_index)
601 consensus_spec[alignment[align_index].first].setIntensity(consensus_spec[alignment[align_index].first].getIntensity() +
602 (pit->getIntensity() / (
double)counter));
604 if (align_index == alignment.size())
609 align_size = align_size + 1 - counter;
614 if (spec_a + spec_b - align_size != consensus_spec.size())
616 LOG_WARN <<
"wrong number of features after merge. Expected: " << spec_a + spec_b - align_size <<
" got: " << consensus_spec.size() <<
"\n";
619 rt_average /= it->second.size() + 1;
620 consensus_spec.
setRT(rt_average);
626 precursor_mz_average /= precursor_count;
631 pcs[0].setMZ(precursor_mz_average);
635 if (consensus_spec.empty())
641 merged_spectra.addSpectrum(consensus_spec);
648 LOG_INFO <<
" size " << it->first <<
": " << it->second <<
"x\n";
652 sprintf(buffer,
"%d/%d (%.2f %%) of blocked spectra", (
int)count_peaks_aligned,
653 (
int)count_peaks_overall,
float(count_peaks_aligned) /
float(count_peaks_overall) * 100.);
659 for (
Size i = 0; i < exp.
size(); ++i)
661 if (merged_indices.count(i) == 0)
676 exp.
getSpectra().insert(exp.
end(), merged_spectra.begin(), merged_spectra.end());
700 template <
typename MapType>
705 double mz_binning_width(param_.getValue(
"mz_binning_width"));
706 String mz_binning_unit(param_.getValue(
"mz_binning_width_unit"));
708 unsigned progress = 0;
709 std::stringstream progress_message;
710 progress_message <<
"averaging profile spectra of MS level " << ms_level;
711 startProgress(0, spectra_to_average_over.size(), progress_message.str());
716 setProgress(++progress);
719 std::vector<double> mz_positions_all;
720 for (std::vector<std::pair<Size, double> >::const_iterator it2 = it->second.begin(); it2 != it->second.end(); ++it2)
725 mz_positions_all.push_back(it_mz->getMZ());
729 sort(mz_positions_all.begin(), mz_positions_all.end());
731 std::vector<double> mz_positions;
732 std::vector<double> intensities;
733 double last_mz = std::numeric_limits<double>::min();
734 double delta_mz(mz_binning_width);
735 for (std::vector<double>::iterator it_mz = mz_positions_all.begin(); it_mz < mz_positions_all.end(); ++it_mz)
737 if (mz_binning_unit ==
"ppm")
739 delta_mz = mz_binning_width * (*it_mz) / 1000000;
742 if (((*it_mz) - last_mz) > delta_mz)
744 mz_positions.push_back(*it_mz);
745 intensities.push_back(0.0);
751 for (std::vector<std::pair<Size, double> >::const_iterator it2 = it->second.begin(); it2 != it->second.end(); ++it2)
757 for (
Size i = 0; i < mz_positions.size(); ++i)
759 if ((spline.
getMzMin() < mz_positions[i]) && (mz_positions[i] < spline.
getMzMax()))
761 intensities[i] += nav.
eval(mz_positions[i]) * (it2->second);
768 average_spec.
clear(
false);
772 for (
Size i = 0; i < mz_positions.size(); ++i)
775 peak.
setMZ(mz_positions[i]);
777 average_spec.push_back(peak);
791 exp[it->first] = exp_tmp[n];
813 template <
typename MapType>
818 double mz_binning_width(param_.getValue(
"mz_binning_width"));
819 String mz_binning_unit(param_.getValue(
"mz_binning_width_unit"));
821 unsigned progress = 0;
823 std::stringstream progress_message;
824 progress_message <<
"averaging centroid spectra of MS level " << ms_level;
825 logger.
startProgress(0, spectra_to_average_over.size(), progress_message.str());
834 std::vector<std::pair<double, double> > mz_intensity_all;
835 for (std::vector<std::pair<Size, double> >::const_iterator it2 = it->second.begin(); it2 != it->second.end(); ++it2)
840 std::pair<double, double> mz_intensity(it_mz->getMZ(), (it_mz->getIntensity() * it2->second));
841 mz_intensity_all.push_back(mz_intensity);
848 std::vector<double> mz_new;
849 std::vector<double> intensity_new;
850 double last_mz = std::numeric_limits<double>::min();
851 double delta_mz = mz_binning_width;
853 double sum_intensity(0);
855 for (std::vector<std::pair<double, double> >::const_iterator it_mz = mz_intensity_all.begin(); it_mz != mz_intensity_all.end(); ++it_mz)
857 if (mz_binning_unit ==
"ppm")
859 delta_mz = mz_binning_width * (it_mz->first) / 1000000;
862 if (((it_mz->first - last_mz) > delta_mz) && (count > 0))
864 mz_new.push_back(sum_mz / count);
865 intensity_new.push_back(sum_intensity);
870 last_mz = it_mz->first;
874 sum_mz += it_mz->first;
875 sum_intensity += it_mz->second;
880 mz_new.push_back(sum_mz / count);
881 intensity_new.push_back(sum_intensity);
886 average_spec.
clear(
false);
890 for (
Size i = 0; i < mz_new.size(); ++i)
893 peak.
setMZ(mz_new[i]);
895 average_spec.push_back(peak);
909 exp[it->first] = exp_tmp[n];
918 bool static compareByFirst(std::pair<double, double> i, std::pair<double, double> j)
920 return i.first < j.first;