From a1c789b89cf984ab64bd52fe1b1017145df0858b Mon Sep 17 00:00:00 2001 From: Joshua Charkow Date: Thu, 17 Apr 2025 17:00:39 -0400 Subject: [PATCH 1/7] first steps of auto IRT refinement: add way to subsample the library --- src/topp/OpenSwathWorkflow.cpp | 55 ++++++++++++++++++++++++++++++++++ 1 file changed, 55 insertions(+) diff --git a/src/topp/OpenSwathWorkflow.cpp b/src/topp/OpenSwathWorkflow.cpp index b183f7a3b3c..5c7cb256b17 100644 --- a/src/topp/OpenSwathWorkflow.cpp +++ b/src/topp/OpenSwathWorkflow.cpp @@ -422,6 +422,9 @@ class TOPPOpenSwathWorkflow registerStringOption_("enable_ms1", "", "true", "Extract the precursor ion trace(s) and use for scoring if present", false, true); setValidStrings_("enable_ms1", ListUtils::create("true,false")); + registerStringOption_("auto_transform", "", "true", "Perform automatic RT, IM, and m/z calibration using a quality heuristic apporach", false, true); + setValidStrings_("auto_transform", ListUtils::create("true,false")); + registerStringOption_("enable_ipf", "", "true", "Enable additional scoring of identification assays using IPF (see online documentation)", false, true); setValidStrings_("enable_ipf", ListUtils::create("true,false")); @@ -663,6 +666,7 @@ class TOPPOpenSwathWorkflow bool pasef = getFlag_("pasef"); bool sort_swath_maps = getFlag_("sort_swath_maps"); bool use_ms1_traces = getStringOption_("enable_ms1") == "true"; + bool auto_transform = getStringOption_("auto_transform") == "true"; bool enable_uis_scoring = getStringOption_("enable_ipf") == "true"; int batchSize = (int)getIntOption_("batchSize"); int outer_loop_threads = (int)getIntOption_("outer_loop_threads"); @@ -904,6 +908,57 @@ class TOPPOpenSwathWorkflow calibration_param.setValue("im_extraction_window", cp_irt.im_extraction_window); calibration_param.setValue("mz_correction_function", mz_correction_function); TransformationDescription trafo_rtnorm; + if (auto_transform) + { + OPENMS_LOG_DEBUG << "Performing automatic RT, IM and m/z calibration" << std::endl; + + // Random subsample of the library to get a smaller set of peptides + auto shuffledCompounds = transition_exp.getCompounds(); + Math::RandomShuffler shuffler_(1); + shuffler_.portable_random_shuffle(shuffledCompounds.begin(), shuffledCompounds.end()); + + // take the first 1000 peptides + OpenSwath::LightTargetedExperiment transition_exp_subsampled; + + // copied from void OpenSwathWorkflow::selectCompoundsForBatch_(const OpenSwath::LightTargetedExperiment& transition_exp_used_all, + // creates a new LightTargetedExperiment from the subsample + // compute batch start/end + int batch_size = 10; + + size_t end = batch_size; + if (end > transition_exp.compounds.size()) + { + end = transition_exp.compounds.size(); + } + + // Create the new, batch-size transition experiment + transition_exp_subsampled.proteins = transition_exp.proteins; + transition_exp_subsampled.compounds.insert(transition_exp_subsampled.compounds.end(), shuffledCompounds.begin(), shuffledCompounds.begin() + end); + + std::set selected_compounds; + for (Size i = 0; i < transition_exp_subsampled.compounds.size(); i++) + { + selected_compounds.insert(transition_exp_subsampled.compounds[i].id); + } + + for (Size i = 0; i < transition_exp.transitions.size(); i++) + { + if (selected_compounds.find(transition_exp.transitions[i].peptide_ref) != selected_compounds.end()) + { + transition_exp_subsampled.transitions.push_back(transition_exp.transitions[i]); + } + } + for (Size i = 0; i < transition_exp_subsampled.getCompounds().size(); i++) + { + OPENMS_LOG_DEBUG << "Subsampled compound " << transition_exp_subsampled.getCompounds()[i].id << std::endl; + } + + for (Size i = 0; i < transition_exp_subsampled.getTransitions().size(); i++) + { + OPENMS_LOG_DEBUG << "Subsampled transition " << transition_exp_subsampled.getTransitions()[i].getNativeID() << std::endl; + } + } + if (nonlinear_irt_tr_file.empty()) { trafo_rtnorm = performCalibration(trafo_in, irt_tr_file, swath_maps, From ccb12f549d50b3ba7bacca109f9174fd162bcfd3 Mon Sep 17 00:00:00 2001 From: Joshua Charkow Date: Fri, 18 Apr 2025 13:10:02 -0400 Subject: [PATCH 2/7] feature: add automatic linear calibration this is not tested yet. Also compute the standard deviation at the end from the predicted amount which will be used as the extraction window for the non linear calibration. --- .../ANALYSIS/OPENSWATH/OpenSwathWorkflow.cpp | 36 +++++++++++++++++-- src/topp/OpenSwathWorkflow.cpp | 28 +++++++++++++++ 2 files changed, 62 insertions(+), 2 deletions(-) diff --git a/src/openms/source/ANALYSIS/OPENSWATH/OpenSwathWorkflow.cpp b/src/openms/source/ANALYSIS/OPENSWATH/OpenSwathWorkflow.cpp index 3693854269b..02f3996d54f 100644 --- a/src/openms/source/ANALYSIS/OPENSWATH/OpenSwathWorkflow.cpp +++ b/src/openms/source/ANALYSIS/OPENSWATH/OpenSwathWorkflow.cpp @@ -209,8 +209,8 @@ namespace OpenMS OPENMS_LOG_DEBUG << "Performed outlier detection, left with features: " << pairs_corrected.size() << std::endl; // 6. Check whether the found peptides fulfill the binned coverage criteria - // set by the user. - if (estimateBestPeptides) + // set by the user. This is only done if not using a linear RT correction + if (estimateBestPeptides && irt_detection_param.getValue("alignmentMethod").toString() != "linear") { bool enoughPeptides = MRMRTNormalizer::computeBinnedCoverage(RTRange, pairs_corrected, irt_detection_param.getValue("NrRTBins"), @@ -274,8 +274,40 @@ namespace OpenMS { OPENMS_LOG_DEBUG << pairs_corrected[i].first << " " << pairs_corrected[i].second << std::endl; } + + /////////////////////////////////////////////////////////////////////////// + // TODO incorperate this into transformation description function instead + ///////////////////////////////////////////////////////////// + + std::vector predicted_values(pairs_corrected.size()); // these are the predicted values by applying the trafo on the values + std::vector delta_true_predicted(pairs_corrected.size()); // these are the differences between the predicted and the true values + double sum_rt_differences = 0.0; + for (size_t i = 0; i < pairs_corrected.size(); ++i) + { + predicted_values[i] = trafo_out.apply(pairs_corrected[i].first); + delta_true_predicted[i] = pairs_corrected[i].second - predicted_values[i]; + sum_rt_differences += delta_true_predicted[i]; // sum for computing the mean + OPENMS_LOG_DEBUG << "True - Predicted value: " << pairs_corrected[i].second << " " << predicted_values[i] << std::endl; + } + double mean_rt_difference = sum_rt_differences / delta_true_predicted.size(); + + // compute the variance + double rt_variance = 0.0; + for (double delta: delta_true_predicted) + { + rt_variance += (delta - mean_rt_difference) * (delta - mean_rt_difference); + } + rt_variance /= delta_true_predicted.size(); + double rt_stdev = sqrt(rt_variance); + double rt_extraction_window = rt_stdev * 3.0 * 2; // 3 times the standard deviation, doubled since the window is the full length (not half) + OPENMS_LOG_DEBUG << "RT variance: " << rt_stdev << std::endl; + OPENMS_LOG_DEBUG << "RT Extraction Window to Use: " << rt_extraction_window << std::endl; + + OPENMS_LOG_DEBUG << "End of doDataNormalization_ method" << std::endl; + trafo_out.printSummary(std::cout); + this->endProgress(); return trafo_out; } diff --git a/src/topp/OpenSwathWorkflow.cpp b/src/topp/OpenSwathWorkflow.cpp index 5c7cb256b17..632c590aca2 100644 --- a/src/topp/OpenSwathWorkflow.cpp +++ b/src/topp/OpenSwathWorkflow.cpp @@ -957,6 +957,34 @@ class TOPPOpenSwathWorkflow { OPENMS_LOG_DEBUG << "Subsampled transition " << transition_exp_subsampled.getTransitions()[i].getNativeID() << std::endl; } + + //////////////////////////////////// + // Do A simple linear calibration from the subsampled library + //////////////////////////////////////////// + std::vector< OpenMS::MSChromatogram > chromatograms; + OpenSwathCalibrationWorkflow wf; + wf.setLogType(log_type_); + wf.simpleExtractChromatograms_(swath_maps, transition_exp_subsampled, chromatograms, + trafo_rtnorm, cp_irt, pasef, load_into_memory); + + // always use estimateBestPeptides for the nonlinear approach + Param irt_lin = irt_detection_param; + + irt_lin.setValue("alignmentMethod", "linear" ); + irt_lin.setValue("estimateBestPeptides", "true"); + + TransformationDescription im_trafo; // exp -> theoretical + trafo_rtnorm = wf.doDataNormalization_(transition_exp_subsampled, chromatograms, im_trafo, swath_maps, + min_rsq, min_coverage, + feature_finder_param, irt_lin, calibration_param, pasef); + + /////////////////////////////////////////////////////////// + // Do a non linear calibration from the subsampled library + /////////////////////////////////////////////////////////// + + Param irt_non_lin = irt_detection_param; + irt_lin.setValue("estimateBestPeptides", "true"); + } if (nonlinear_irt_tr_file.empty()) From 300e9998389a603f2ad04d4082d6337d20357c53 Mon Sep 17 00:00:00 2001 From: Joshua Charkow Date: Wed, 30 Apr 2025 08:56:12 -0400 Subject: [PATCH 3/7] add method to subsample LightTargetedExperiment based on peptide list --- .../DATAACCESS/TransitionExperiment.h | 24 +++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/src/openswathalgo/include/OpenMS/OPENSWATHALGO/DATAACCESS/TransitionExperiment.h b/src/openswathalgo/include/OpenMS/OPENSWATHALGO/DATAACCESS/TransitionExperiment.h index 9f6ebe5c4bb..78a4c525ac2 100644 --- a/src/openswathalgo/include/OpenMS/OPENSWATHALGO/DATAACCESS/TransitionExperiment.h +++ b/src/openswathalgo/include/OpenMS/OPENSWATHALGO/DATAACCESS/TransitionExperiment.h @@ -238,6 +238,30 @@ namespace OpenSwath return *(compound_reference_map_[ref]); } + // Creates a subsampled library from the specified light compounds + const LightTargetedExperiment selectCompounds(const std::vector& selectedCompounds) + { + LightTargetedExperiment subsampled; + subsampled.proteins = proteins; + + subsampled.compounds.insert(subsampled.compounds.end(), selectedCompounds.begin(), selectedCompounds.end()); + + std::set compound_ids; + for (size_t i = 0; i < subsampled.compounds.size(); i++) + { + compound_ids.insert(subsampled.compounds[i].id); + } + + for (size_t i = 0; i < transitions.size(); i++) + { + if (compound_ids.find(transitions[i].peptide_ref) != compound_ids.end()) + { + subsampled.transitions.push_back(subsampled.transitions[i]); + } + } + return subsampled; + } + private: void createPeptideReferenceMap_() From 4d5fa96ddadf01a7dc9d695a49596eb050e90b58 Mon Sep 17 00:00:00 2001 From: Joshua Charkow Date: Wed, 30 Apr 2025 09:30:31 -0400 Subject: [PATCH 4/7] fix, refactor: subsampleLibrary across RT space --- .../ANALYSIS/OPENSWATH/OpenSwathHelper.h | 16 ++++++ .../ANALYSIS/OPENSWATH/OpenSwathHelper.cpp | 55 +++++++++++++++++++ .../DATAACCESS/TransitionExperiment.h | 3 +- 3 files changed, 73 insertions(+), 1 deletion(-) diff --git a/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/OpenSwathHelper.h b/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/OpenSwathHelper.h index 4cc25f44ec0..d4d78b15df3 100644 --- a/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/OpenSwathHelper.h +++ b/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/OpenSwathHelper.h @@ -198,6 +198,22 @@ namespace OpenMS static std::map simpleFindBestFeature(const OpenMS::MRMFeatureFinderScoring::TransitionGroupMapType & transition_group_map, bool useQualCutoff = false, double qualCutoff = 0.0); + + /** + * @brief Subsample a library based on a given + * + * This method is used to subsample the library intelligently in order to ensure that the subsampled compounds are distributed across RT space. + * + * @param nrBins Number of bins across RT + * @param peptidesPerBin aim number of peptides per bin, if not reached a warning is issued + * + * @return LightTargetedExperiment with the subsampled library + */ + OpenSwath::LightTargetedExperiment subsampleLibrary(OpenSwath::LightTargetedExperiment& library, + int nrBins, + int minPeptidesPerBin); + + }; } // namespace OpenMS diff --git a/src/openms/source/ANALYSIS/OPENSWATH/OpenSwathHelper.cpp b/src/openms/source/ANALYSIS/OPENSWATH/OpenSwathHelper.cpp index f810215e2c7..4a2b8229760 100644 --- a/src/openms/source/ANALYSIS/OPENSWATH/OpenSwathHelper.cpp +++ b/src/openms/source/ANALYSIS/OPENSWATH/OpenSwathHelper.cpp @@ -7,6 +7,7 @@ // -------------------------------------------------------------------------- #include +#include namespace OpenMS { @@ -187,4 +188,58 @@ namespace OpenMS return result; } + OpenSwath::LightTargetedExperiment OpenSwathHelper::subsampleLibrary(OpenSwath::LightTargetedExperiment& library, + int nrBins, + int peptidesPerBin) + { + std::pair RTRange = OpenSwathHelper::estimateRTRange(library); + OPENMS_LOG_DEBUG << "Detected retention time range from " << RTRange.first << " to " << RTRange.second << std::endl; + + // Sort the LightTargetedExperiment Peptide precursors by their RT + auto peptide_precursors = library.getCompounds(); + std::sort(peptide_precursors.begin(), peptide_precursors.end(), + [](const OpenSwath::LightCompound& a, const OpenSwath::LightCompound& b) + { + return a.rt < b.rt; + }); + + std::vector selected_compounds; + + // Create a histogram bins with the desired number of bins + Math::Histogram<> hist(RTRange.first, RTRange.second, nrBins); + + // Take the first `minPeptidesPerBin` peptides from each bin + for (size_t h = 0; h < hist.size(); h++) + { + int numPepsInBin = 0; + int i = 0; + bool doneWithBin = false; + while ((numPepsInBin < peptidesPerBin) && (!doneWithBin)) + { + // Check if the peptide precursor is within the bin range + if (peptide_precursors[i].rt >= hist.leftBorderOfBin(h) && peptide_precursors[i].rt < hist.rightBorderOfBin(h)) + { + numPepsInBin++; + i++; + selected_compounds.push_back(peptide_precursors[i]); + } + else + { + OPENMS_LOG_WARN << "RT Bin from" << hist.leftBorderOfBin(h) << " to " << hist.rightBorderOfBin(h + 1) << " has only" + << numPepsInBin<< " peptides. This is less than the minumum number specified of " << peptidesPerBin << std::endl; + + // Have the minumum number of peptides required, seek to the next bin + while (peptide_precursors[i].rt < hist.rightBorderOfBin(h)) + { + i++; + } + doneWithBin = true; + numPepsInBin = 0; + } + } + } + return library.selectCompounds(selected_compounds); + } + + } diff --git a/src/openswathalgo/include/OpenMS/OPENSWATHALGO/DATAACCESS/TransitionExperiment.h b/src/openswathalgo/include/OpenMS/OPENSWATHALGO/DATAACCESS/TransitionExperiment.h index 78a4c525ac2..682e503f954 100644 --- a/src/openswathalgo/include/OpenMS/OPENSWATHALGO/DATAACCESS/TransitionExperiment.h +++ b/src/openswathalgo/include/OpenMS/OPENSWATHALGO/DATAACCESS/TransitionExperiment.h @@ -11,6 +11,7 @@ #include #include #include +#include #include @@ -239,7 +240,7 @@ namespace OpenSwath } // Creates a subsampled library from the specified light compounds - const LightTargetedExperiment selectCompounds(const std::vector& selectedCompounds) + LightTargetedExperiment selectCompounds(const std::vector& selectedCompounds) { LightTargetedExperiment subsampled; subsampled.proteins = proteins; From 5526f70564b9575f2d6ffb8bf6a5007ee968c8e3 Mon Sep 17 00:00:00 2001 From: Joshua Charkow Date: Wed, 30 Apr 2025 10:40:15 -0400 Subject: [PATCH 5/7] initial framework for auto irt detection with linear transformation --- .../ANALYSIS/OPENSWATH/OpenSwathHelper.h | 2 +- .../ANALYSIS/OPENSWATH/OpenSwathWorkflow.h | 15 ++ .../ANALYSIS/OPENSWATH/OpenSwathWorkflow.cpp | 210 ++++++++++++++++++ src/topp/OpenSwathWorkflow.cpp | 2 +- 4 files changed, 227 insertions(+), 2 deletions(-) diff --git a/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/OpenSwathHelper.h b/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/OpenSwathHelper.h index d4d78b15df3..4496559b656 100644 --- a/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/OpenSwathHelper.h +++ b/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/OpenSwathHelper.h @@ -209,7 +209,7 @@ namespace OpenMS * * @return LightTargetedExperiment with the subsampled library */ - OpenSwath::LightTargetedExperiment subsampleLibrary(OpenSwath::LightTargetedExperiment& library, + static OpenSwath::LightTargetedExperiment subsampleLibrary(OpenSwath::LightTargetedExperiment& library, int nrBins, int minPeptidesPerBin); diff --git a/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/OpenSwathWorkflow.h b/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/OpenSwathWorkflow.h index a6622d72e7d..5c09daa3712 100644 --- a/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/OpenSwathWorkflow.h +++ b/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/OpenSwathWorkflow.h @@ -280,6 +280,21 @@ namespace OpenMS bool pasef = false, bool load_into_memory = false); + TransformationDescription performAutoRTNormalization( + OpenSwath::LightTargetedExperiment& irt_transitions, + std::vector< OpenSwath::SwathMap > & swath_maps, + TransformationDescription& im_trafo, + double min_rsq, + double min_coverage, + const Param& feature_finder_param, + const ChromExtractParams& cp_irt, + const Param& irt_detection_param, + const Param& calibration_param, + const String& irt_mzml_out, + Size debug_level, + bool pasef, + bool load_into_memory); + public: /** @brief Perform retention time and m/z calibration diff --git a/src/openms/source/ANALYSIS/OPENSWATH/OpenSwathWorkflow.cpp b/src/openms/source/ANALYSIS/OPENSWATH/OpenSwathWorkflow.cpp index 02f3996d54f..bf4457f5bb7 100644 --- a/src/openms/source/ANALYSIS/OPENSWATH/OpenSwathWorkflow.cpp +++ b/src/openms/source/ANALYSIS/OPENSWATH/OpenSwathWorkflow.cpp @@ -86,6 +86,216 @@ namespace OpenMS return tr; } + TransformationDescription OpenSwathCalibrationWorkflow::performAutoRTNormalization( + OpenSwath::LightTargetedExperiment& targeted_exp, + std::vector< OpenSwath::SwathMap > & swath_maps, + TransformationDescription& im_trafo, + double min_rsq, + double min_coverage, + const Param& default_ffparam, + const ChromExtractParams& cp_irt, + const Param& irt_detection_param, + const Param& calibration_param, + const String& irt_mzml_out, + Size debug_level, + bool pasef, + bool load_into_memory) + { + OPENMS_LOG_DEBUG << "Start of performAutoRTNormalization method" << std::endl; + + // 1. Estimate the retention time range of the iRT peptides over all assays + std::pair RTRange = OpenSwathHelper::estimateRTRange(targeted_exp); + OPENMS_LOG_DEBUG << "Detected retention time range from " << RTRange.first << " to " << RTRange.second << std::endl; + + // ########## LINEAR CALIBRATION ########## + // 1. Subsample the library to a number of bins + OpenSwath::LightTargetedExperiment targeted_exp_subsampled = OpenSwathHelper::subsampleLibrary(targeted_exp, + irt_detection_param.getValue("numberOfBins"), + irt_detection_param.getValue("peptidesPerBin")); + + // 2. Store the peptide retention times in an intermediate map + std::map PeptideRTMap; + for (Size i = 0; i < targeted_exp_subsampled.getCompounds().size(); i++) + { + PeptideRTMap[targeted_exp_subsampled.getCompounds()[i].id] = targeted_exp_subsampled.getCompounds()[i].rt; + } + + // 2. Extract the chromatograms + std::vector< OpenMS::MSChromatogram > chromatograms; + TransformationDescription trafo; // dummy + this->simpleExtractChromatograms_(swath_maps, targeted_exp_subsampled, chromatograms, trafo, cp_irt, pasef, load_into_memory); + + OPENMS_LOG_DEBUG << "Extracted number of chromatograms from iRT files: " << chromatograms.size() << std::endl; + + // 2. For each peptide, extract and score chromatograms + const OpenSwath::LightTargetedExperiment& transition_exp_used = targeted_exp_subsampled; + + // Change the feature finding parameters: + // - no RT score (since we don't know the correct retention time) + // - no RT window + // - no elution model score + // - no peak quality (use all peaks) + // - if best peptides should be used, use peak quality + MRMFeatureFinderScoring featureFinder; + Param feature_finder_param(default_ffparam); + feature_finder_param.setValue("Scores:use_rt_score", "false"); + feature_finder_param.setValue("Scores:use_elution_model_score", "false"); + feature_finder_param.setValue("rt_extraction_window", -1.0); + feature_finder_param.setValue("stop_report_after_feature", 1); + feature_finder_param.setValue("TransitionGroupPicker:PeakPickerChromatogram:signal_to_noise", 1.0); // set to 1.0 in all cases + feature_finder_param.setValue("TransitionGroupPicker:compute_peak_quality", "false"); // no peak quality -> take all peaks! + feature_finder_param.setValue("TransitionGroupPicker:compute_peak_quality", "true"); + feature_finder_param.setValue("TransitionGroupPicker:minimal_quality", irt_detection_param.getValue("InitialQualityCutoff")); + featureFinder.setParameters(feature_finder_param); + + FeatureMap featureFile; // for results + OpenMS::MRMFeatureFinderScoring::TransitionGroupMapType transition_group_map; // for results + std::vector empty_swath_maps; + TransformationDescription empty_trafo; // empty transformation + + // Prepare the data with the chromatograms + boost::shared_ptr xic_map(new PeakMap); + xic_map->setChromatograms(chromatograms); + OpenSwath::SpectrumAccessPtr chromatogram_ptr = OpenSwath::SpectrumAccessPtr(new OpenMS::SpectrumAccessOpenMS(xic_map)); + + featureFinder.pickExperiment(chromatogram_ptr, featureFile, transition_exp_used, empty_trafo, empty_swath_maps, transition_group_map); + + // 4. Find most likely correct feature for each compound and add it to the + // "pairs" vector by computing pairs of iRT and real RT. + // + // Note that the quality threshold will only be applied if + // estimateBestPeptides is true + std::vector > pairs; // store the RT pairs to write the output trafoXML + std::map best_features = OpenSwathHelper::simpleFindBestFeature(transition_group_map, + true, irt_detection_param.getValue("OverallQualityCutoff")); + OPENMS_LOG_DEBUG << "Extracted best features: " << best_features.size() << std::endl; + + // Create pairs vector and store peaks + std::map trgrmap_allpeaks; // store all peaks above cutoff + for (std::map::iterator it = best_features.begin(); it != best_features.end(); ++it) + { + pairs.emplace_back(it->second, PeptideRTMap[it->first]); // pair + if (transition_group_map.find(it->first) != transition_group_map.end()) + { + trgrmap_allpeaks[ it->first ] = &transition_group_map[ it->first]; + } + } + + // 5. Perform the outlier detection + std::vector > pairs_corrected; + String outlier_method = irt_detection_param.getValue("outlierMethod").toString(); + if (outlier_method == "iter_residual" || outlier_method == "iter_jackknife") + { + pairs_corrected = MRMRTNormalizer::removeOutliersIterative(pairs, min_rsq, min_coverage, + irt_detection_param.getValue("useIterativeChauvenet").toBool(), outlier_method); + } + else if (outlier_method == "ransac") + { + // First, estimate of the maximum deviation from RT that is tolerated: + // Because 120 min gradient can have around 4 min elution shift, we use + // a default value of 3 % of the gradient to find upper RT threshold (3.6 min). + double pcnt_rt_threshold = irt_detection_param.getValue("RANSACMaxPercentRTThreshold"); + double max_rt_threshold = (RTRange.second - RTRange.first) * pcnt_rt_threshold / 100.0; + + pairs_corrected = MRMRTNormalizer::removeOutliersRANSAC(pairs, min_rsq, min_coverage, + irt_detection_param.getValue("RANSACMaxIterations"), max_rt_threshold, + irt_detection_param.getValue("RANSACSamplingSize")); + } + else if (outlier_method == "none") + { + pairs_corrected = pairs; + } + else + { + throw Exception::IllegalArgument(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, + String("Illegal argument '") + outlier_method + + "' used for outlierMethod (valid: 'iter_residual', 'iter_jackknife', 'ransac', 'none')."); + } + OPENMS_LOG_DEBUG << "Performed outlier detection, left with features: " << pairs_corrected.size() << std::endl; + + // 7. Select the "correct" peaks for m/z (and IM) correction (e.g. remove those not + // part of the linear regression) + std::map trgrmap_final; // store all peaks above cutoff + for (const auto& it : trgrmap_allpeaks) + { + if (it.second->getFeatures().empty() ) {continue;} + const MRMFeature& feat = it.second->getBestFeature(); + + // Check if the current feature is in the list of pairs used for the + // linear RT regression (using other features may result in wrong + // calibration values). + // Matching only by RT is not perfect but should work for most cases. + for (Size pit = 0; pit < pairs_corrected.size(); pit++) + { + if (fabs(feat.getRT() - pairs_corrected[pit].first ) < 1e-2) + { + trgrmap_final[ it.first ] = it.second; + break; + } + } + } + + // 8. Correct m/z (and IM) deviations using SwathMapMassCorrection + // m/z correction is done with the -irt_im_extraction parameters + SwathMapMassCorrection mc; + mc.setParameters(calibration_param); + + mc.correctMZ(trgrmap_final, targeted_exp, swath_maps, pasef); + mc.correctIM(trgrmap_final, targeted_exp, swath_maps, pasef, im_trafo); + + // 9. store RT transformation, using the selected model + TransformationDescription trafo_out; + trafo_out.setDataPoints(pairs_corrected); + Param model_params; + model_params.setValue("symmetric_regression", "false"); + model_params.setValue("span", irt_detection_param.getValue("lowess:span")); + model_params.setValue("num_nodes", irt_detection_param.getValue("b_spline:num_nodes")); + String model_type = irt_detection_param.getValue("alignmentMethod").toString(); + trafo_out.fitModel(model_type, model_params); + + OPENMS_LOG_DEBUG << "Final RT mapping:" << std::endl; + for (Size i = 0; i < pairs_corrected.size(); i++) + { + OPENMS_LOG_DEBUG << pairs_corrected[i].first << " " << pairs_corrected[i].second << std::endl; + } + + /////////////////////////////////////////////////////////////////////////// + // TODO incorperate this into transformation description function instead + ///////////////////////////////////////////////////////////// + + std::vector predicted_values(pairs_corrected.size()); // these are the predicted values by applying the trafo on the values + std::vector delta_true_predicted(pairs_corrected.size()); // these are the differences between the predicted and the true values + double sum_rt_differences = 0.0; + for (size_t i = 0; i < pairs_corrected.size(); ++i) + { + predicted_values[i] = trafo_out.apply(pairs_corrected[i].first); + delta_true_predicted[i] = pairs_corrected[i].second - predicted_values[i]; + sum_rt_differences += delta_true_predicted[i]; // sum for computing the mean + OPENMS_LOG_DEBUG << "True - Predicted value: " << pairs_corrected[i].second << " " << predicted_values[i] << std::endl; + } + double mean_rt_difference = sum_rt_differences / delta_true_predicted.size(); + + // compute the variance + double rt_variance = 0.0; + for (double delta: delta_true_predicted) + { + rt_variance += (delta - mean_rt_difference) * (delta - mean_rt_difference); + } + rt_variance /= delta_true_predicted.size(); + double rt_stdev = sqrt(rt_variance); + double rt_extraction_window = rt_stdev * 3.0 * 2; // 3 times the standard deviation, doubled since the window is the full length (not half) + OPENMS_LOG_DEBUG << "RT variance: " << rt_stdev << std::endl; + OPENMS_LOG_DEBUG << "RT Extraction Window to Use: " << rt_extraction_window << std::endl; + + + OPENMS_LOG_DEBUG << "End of doDataNormalization_ method" << std::endl; + + trafo_out.printSummary(std::cout); + + this->endProgress(); + return trafo_out; + } + TransformationDescription OpenSwathCalibrationWorkflow::doDataNormalization_( const OpenSwath::LightTargetedExperiment& targeted_exp, const std::vector< OpenMS::MSChromatogram >& chromatograms, diff --git a/src/topp/OpenSwathWorkflow.cpp b/src/topp/OpenSwathWorkflow.cpp index 632c590aca2..5a3298bf5f8 100644 --- a/src/topp/OpenSwathWorkflow.cpp +++ b/src/topp/OpenSwathWorkflow.cpp @@ -967,7 +967,7 @@ class TOPPOpenSwathWorkflow wf.simpleExtractChromatograms_(swath_maps, transition_exp_subsampled, chromatograms, trafo_rtnorm, cp_irt, pasef, load_into_memory); - // always use estimateBestPeptides for the nonlinear approach + Param irt_lin = irt_detection_param; irt_lin.setValue("alignmentMethod", "linear" ); From 9e369d910de07364e7cf7826b306e0c23adf7aab Mon Sep 17 00:00:00 2001 From: Joshua Charkow Date: Wed, 30 Apr 2025 17:16:50 -0400 Subject: [PATCH 6/7] progress: auto iRT fix runtime bugs. have a working implementation of the linear auto irt refinement! --- .../ANALYSIS/OPENSWATH/OpenSwathWorkflow.h | 1 - .../ANALYSIS/OPENSWATH/OpenSwathHelper.cpp | 30 ++++--- .../ANALYSIS/OPENSWATH/OpenSwathWorkflow.cpp | 7 +- .../DATAACCESS/TransitionExperiment.h | 4 +- src/topp/OpenSwathWorkflow.cpp | 81 +++---------------- 5 files changed, 36 insertions(+), 87 deletions(-) diff --git a/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/OpenSwathWorkflow.h b/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/OpenSwathWorkflow.h index 5c09daa3712..be402f38c80 100644 --- a/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/OpenSwathWorkflow.h +++ b/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/OpenSwathWorkflow.h @@ -290,7 +290,6 @@ namespace OpenMS const ChromExtractParams& cp_irt, const Param& irt_detection_param, const Param& calibration_param, - const String& irt_mzml_out, Size debug_level, bool pasef, bool load_into_memory); diff --git a/src/openms/source/ANALYSIS/OPENSWATH/OpenSwathHelper.cpp b/src/openms/source/ANALYSIS/OPENSWATH/OpenSwathHelper.cpp index 4a2b8229760..82c76c633d5 100644 --- a/src/openms/source/ANALYSIS/OPENSWATH/OpenSwathHelper.cpp +++ b/src/openms/source/ANALYSIS/OPENSWATH/OpenSwathHelper.cpp @@ -207,37 +207,47 @@ namespace OpenMS // Create a histogram bins with the desired number of bins Math::Histogram<> hist(RTRange.first, RTRange.second, nrBins); + for (size_t h = 0; h < hist.size(); h++) + { + OPENMS_LOG_DEBUG << "Bin " << h << " has left border of " << hist.leftBorderOfBin(h) << " and right border of " << hist.rightBorderOfBin(h) << std::endl; + } // Take the first `minPeptidesPerBin` peptides from each bin + int i = 0; // index of the peptide precursor for (size_t h = 0; h < hist.size(); h++) { int numPepsInBin = 0; - int i = 0; bool doneWithBin = false; while ((numPepsInBin < peptidesPerBin) && (!doneWithBin)) { // Check if the peptide precursor is within the bin range + OPENMS_LOG_DEBUG << "Checking peptide precursor with RT " << peptide_precursors[i].rt << " against bin with left border of " << hist.leftBorderOfBin(h) << std::endl; if (peptide_precursors[i].rt >= hist.leftBorderOfBin(h) && peptide_precursors[i].rt < hist.rightBorderOfBin(h)) { + OPENMS_LOG_DEBUG << "Adding peptide precursor with RT " << peptide_precursors[i].rt << " to bin with left border of " << hist.leftBorderOfBin(h) << std::endl; + selected_compounds.push_back(peptide_precursors[i]); numPepsInBin++; i++; - selected_compounds.push_back(peptide_precursors[i]); } else { - OPENMS_LOG_WARN << "RT Bin from" << hist.leftBorderOfBin(h) << " to " << hist.rightBorderOfBin(h + 1) << " has only" - << numPepsInBin<< " peptides. This is less than the minumum number specified of " << peptidesPerBin << std::endl; + // If we have less than the minimum number of peptides in the bin, warn the user + OPENMS_LOG_WARN << "RT Bin from " << hist.leftBorderOfBin(h) << " to " << hist.rightBorderOfBin(h) << " has only " + << numPepsInBin << " peptides. This is less than the minumum number specified of " << peptidesPerBin << std::endl; - // Have the minumum number of peptides required, seek to the next bin - while (peptide_precursors[i].rt < hist.rightBorderOfBin(h)) - { - i++; - } doneWithBin = true; - numPepsInBin = 0; } } + + // Have the minumum number of peptides required, seek to the next bin + while (peptide_precursors[i].rt < hist.rightBorderOfBin(h)) + { + i++; + } + numPepsInBin = 0; } + + OPENMS_LOG_DEBUG << "Subsampled library has " << selected_compounds.size() << " precursors" << std::endl; return library.selectCompounds(selected_compounds); } diff --git a/src/openms/source/ANALYSIS/OPENSWATH/OpenSwathWorkflow.cpp b/src/openms/source/ANALYSIS/OPENSWATH/OpenSwathWorkflow.cpp index bf4457f5bb7..53ed1a277a5 100644 --- a/src/openms/source/ANALYSIS/OPENSWATH/OpenSwathWorkflow.cpp +++ b/src/openms/source/ANALYSIS/OPENSWATH/OpenSwathWorkflow.cpp @@ -96,7 +96,6 @@ namespace OpenMS const ChromExtractParams& cp_irt, const Param& irt_detection_param, const Param& calibration_param, - const String& irt_mzml_out, Size debug_level, bool pasef, bool load_into_memory) @@ -109,9 +108,10 @@ namespace OpenMS // ########## LINEAR CALIBRATION ########## // 1. Subsample the library to a number of bins + // TODO change param names NrRTBins and MinPeptidesPerBin to more descriptive names for new workflow (e.g. MinPeptidesPerBin is not actually bin peptides, NrRTBins is ok) OpenSwath::LightTargetedExperiment targeted_exp_subsampled = OpenSwathHelper::subsampleLibrary(targeted_exp, - irt_detection_param.getValue("numberOfBins"), - irt_detection_param.getValue("peptidesPerBin")); + irt_detection_param.getValue("NrRTBins"), + irt_detection_param.getValue("MinPeptidesPerBin")); // 2. Store the peptide retention times in an intermediate map std::map PeptideRTMap; @@ -147,6 +147,7 @@ namespace OpenMS feature_finder_param.setValue("TransitionGroupPicker:compute_peak_quality", "true"); feature_finder_param.setValue("TransitionGroupPicker:minimal_quality", irt_detection_param.getValue("InitialQualityCutoff")); featureFinder.setParameters(feature_finder_param); + featureFinder.setStrictFlag(false); // Since we do not know if features are correct, we should not be strict FeatureMap featureFile; // for results OpenMS::MRMFeatureFinderScoring::TransitionGroupMapType transition_group_map; // for results diff --git a/src/openswathalgo/include/OpenMS/OPENSWATHALGO/DATAACCESS/TransitionExperiment.h b/src/openswathalgo/include/OpenMS/OPENSWATHALGO/DATAACCESS/TransitionExperiment.h index 682e503f954..33b34d03a88 100644 --- a/src/openswathalgo/include/OpenMS/OPENSWATHALGO/DATAACCESS/TransitionExperiment.h +++ b/src/openswathalgo/include/OpenMS/OPENSWATHALGO/DATAACCESS/TransitionExperiment.h @@ -240,7 +240,7 @@ namespace OpenSwath } // Creates a subsampled library from the specified light compounds - LightTargetedExperiment selectCompounds(const std::vector& selectedCompounds) + const LightTargetedExperiment selectCompounds(const std::vector& selectedCompounds) { LightTargetedExperiment subsampled; subsampled.proteins = proteins; @@ -257,7 +257,7 @@ namespace OpenSwath { if (compound_ids.find(transitions[i].peptide_ref) != compound_ids.end()) { - subsampled.transitions.push_back(subsampled.transitions[i]); + subsampled.transitions.push_back(transitions[i]); } } return subsampled; diff --git a/src/topp/OpenSwathWorkflow.cpp b/src/topp/OpenSwathWorkflow.cpp index 5a3298bf5f8..1dc90818719 100644 --- a/src/topp/OpenSwathWorkflow.cpp +++ b/src/topp/OpenSwathWorkflow.cpp @@ -910,83 +910,22 @@ class TOPPOpenSwathWorkflow TransformationDescription trafo_rtnorm; if (auto_transform) { - OPENMS_LOG_DEBUG << "Performing automatic RT, IM and m/z calibration" << std::endl; - // Random subsample of the library to get a smaller set of peptides - auto shuffledCompounds = transition_exp.getCompounds(); - Math::RandomShuffler shuffler_(1); - shuffler_.portable_random_shuffle(shuffledCompounds.begin(), shuffledCompounds.end()); - - // take the first 1000 peptides - OpenSwath::LightTargetedExperiment transition_exp_subsampled; - - // copied from void OpenSwathWorkflow::selectCompoundsForBatch_(const OpenSwath::LightTargetedExperiment& transition_exp_used_all, - // creates a new LightTargetedExperiment from the subsample - // compute batch start/end - int batch_size = 10; - - size_t end = batch_size; - if (end > transition_exp.compounds.size()) - { - end = transition_exp.compounds.size(); - } - - // Create the new, batch-size transition experiment - transition_exp_subsampled.proteins = transition_exp.proteins; - transition_exp_subsampled.compounds.insert(transition_exp_subsampled.compounds.end(), shuffledCompounds.begin(), shuffledCompounds.begin() + end); - - std::set selected_compounds; - for (Size i = 0; i < transition_exp_subsampled.compounds.size(); i++) - { - selected_compounds.insert(transition_exp_subsampled.compounds[i].id); - } - - for (Size i = 0; i < transition_exp.transitions.size(); i++) - { - if (selected_compounds.find(transition_exp.transitions[i].peptide_ref) != selected_compounds.end()) - { - transition_exp_subsampled.transitions.push_back(transition_exp.transitions[i]); - } - } - for (Size i = 0; i < transition_exp_subsampled.getCompounds().size(); i++) - { - OPENMS_LOG_DEBUG << "Subsampled compound " << transition_exp_subsampled.getCompounds()[i].id << std::endl; - } - - for (Size i = 0; i < transition_exp_subsampled.getTransitions().size(); i++) - { - OPENMS_LOG_DEBUG << "Subsampled transition " << transition_exp_subsampled.getTransitions()[i].getNativeID() << std::endl; - } - - //////////////////////////////////// - // Do A simple linear calibration from the subsampled library - //////////////////////////////////////////// - std::vector< OpenMS::MSChromatogram > chromatograms; OpenSwathCalibrationWorkflow wf; wf.setLogType(log_type_); - wf.simpleExtractChromatograms_(swath_maps, transition_exp_subsampled, chromatograms, - trafo_rtnorm, cp_irt, pasef, load_into_memory); - - - Param irt_lin = irt_detection_param; + Param linear_irt = irt_detection_param; + linear_irt.setValue("alignmentMethod", "linear"); + Param no_calibration = calibration_param; + no_calibration.setValue("mz_correction_function", "none"); - irt_lin.setValue("alignmentMethod", "linear" ); - irt_lin.setValue("estimateBestPeptides", "true"); + std::vector< OpenMS::MSChromatogram > chromatograms; TransformationDescription im_trafo; // exp -> theoretical - trafo_rtnorm = wf.doDataNormalization_(transition_exp_subsampled, chromatograms, im_trafo, swath_maps, - min_rsq, min_coverage, - feature_finder_param, irt_lin, calibration_param, pasef); - - /////////////////////////////////////////////////////////// - // Do a non linear calibration from the subsampled library - /////////////////////////////////////////////////////////// - - Param irt_non_lin = irt_detection_param; - irt_lin.setValue("estimateBestPeptides", "true"); - - } - + trafo_rtnorm = wf.performAutoRTNormalization(transition_exp, swath_maps, im_trafo, + min_rsq, min_coverage, feature_finder_param, + cp_irt, linear_irt, no_calibration, + debug_level, pasef, load_into_memory); + } if (nonlinear_irt_tr_file.empty()) { trafo_rtnorm = performCalibration(trafo_in, irt_tr_file, swath_maps, From f062df3922af814de30bf7a706f11b8704ba5878 Mon Sep 17 00:00:00 2001 From: Joshua Charkow Date: Fri, 16 May 2025 13:44:12 -0400 Subject: [PATCH 7/7] feature: add auto non linear irt correction --- .../ANALYSIS/OPENSWATH/OpenSwathWorkflow.h | 2 +- .../ANALYSIS/OPENSWATH/OpenSwathWorkflow.cpp | 529 ++++++++++++------ src/topp/OpenSwathWorkflow.cpp | 6 + 3 files changed, 374 insertions(+), 163 deletions(-) diff --git a/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/OpenSwathWorkflow.h b/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/OpenSwathWorkflow.h index be402f38c80..4f738991645 100644 --- a/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/OpenSwathWorkflow.h +++ b/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/OpenSwathWorkflow.h @@ -287,7 +287,7 @@ namespace OpenMS double min_rsq, double min_coverage, const Param& feature_finder_param, - const ChromExtractParams& cp_irt, + ChromExtractParams& cp_irt, const Param& irt_detection_param, const Param& calibration_param, Size debug_level, diff --git a/src/openms/source/ANALYSIS/OPENSWATH/OpenSwathWorkflow.cpp b/src/openms/source/ANALYSIS/OPENSWATH/OpenSwathWorkflow.cpp index 53ed1a277a5..ffeafa5b324 100644 --- a/src/openms/source/ANALYSIS/OPENSWATH/OpenSwathWorkflow.cpp +++ b/src/openms/source/ANALYSIS/OPENSWATH/OpenSwathWorkflow.cpp @@ -93,208 +93,413 @@ namespace OpenMS double min_rsq, double min_coverage, const Param& default_ffparam, - const ChromExtractParams& cp_irt, + ChromExtractParams& cp_irt, const Param& irt_detection_param, const Param& calibration_param, Size debug_level, bool pasef, bool load_into_memory) { - OPENMS_LOG_DEBUG << "Start of performAutoRTNormalization method" << std::endl; + int rt_extraction_window; + TransformationDescription trafo_out; // 1. Estimate the retention time range of the iRT peptides over all assays std::pair RTRange = OpenSwathHelper::estimateRTRange(targeted_exp); OPENMS_LOG_DEBUG << "Detected retention time range from " << RTRange.first << " to " << RTRange.second << std::endl; - // ########## LINEAR CALIBRATION ########## - // 1. Subsample the library to a number of bins - // TODO change param names NrRTBins and MinPeptidesPerBin to more descriptive names for new workflow (e.g. MinPeptidesPerBin is not actually bin peptides, NrRTBins is ok) - OpenSwath::LightTargetedExperiment targeted_exp_subsampled = OpenSwathHelper::subsampleLibrary(targeted_exp, - irt_detection_param.getValue("NrRTBins"), - irt_detection_param.getValue("MinPeptidesPerBin")); - - // 2. Store the peptide retention times in an intermediate map - std::map PeptideRTMap; - for (Size i = 0; i < targeted_exp_subsampled.getCompounds().size(); i++) { - PeptideRTMap[targeted_exp_subsampled.getCompounds()[i].id] = targeted_exp_subsampled.getCompounds()[i].rt; - } - - // 2. Extract the chromatograms - std::vector< OpenMS::MSChromatogram > chromatograms; - TransformationDescription trafo; // dummy - this->simpleExtractChromatograms_(swath_maps, targeted_exp_subsampled, chromatograms, trafo, cp_irt, pasef, load_into_memory); - - OPENMS_LOG_DEBUG << "Extracted number of chromatograms from iRT files: " << chromatograms.size() << std::endl; - - // 2. For each peptide, extract and score chromatograms - const OpenSwath::LightTargetedExperiment& transition_exp_used = targeted_exp_subsampled; + OPENMS_LOG_DEBUG << "Start of performAutoRTNormalization method" << std::endl; + + // ########## LINEAR CALIBRATION ########## + // 1. Subsample the library to a number of bins + // TODO change param names NrRTBins and MinPeptidesPerBin to more descriptive names for new workflow (e.g. MinPeptidesPerBin is not actually bin peptides, NrRTBins is ok) + OpenSwath::LightTargetedExperiment targeted_exp_subsampled = OpenSwathHelper::subsampleLibrary(targeted_exp, + irt_detection_param.getValue("NrRTBins"), + irt_detection_param.getValue("MinPeptidesPerBin")); + + // 2. Store the peptide retention times in an intermediate map + std::map PeptideRTMap; + for (Size i = 0; i < targeted_exp_subsampled.getCompounds().size(); i++) + { + PeptideRTMap[targeted_exp_subsampled.getCompounds()[i].id] = targeted_exp_subsampled.getCompounds()[i].rt; + } + + // 2. Extract the chromatograms + std::vector< OpenMS::MSChromatogram > chromatograms; + TransformationDescription trafo; // dummy + this->simpleExtractChromatograms_(swath_maps, targeted_exp_subsampled, chromatograms, trafo, cp_irt, pasef, load_into_memory); + + OPENMS_LOG_DEBUG << "Extracted number of chromatograms from iRT files: " << chromatograms.size() << std::endl; + + // 2. For each peptide, extract and score chromatograms + const OpenSwath::LightTargetedExperiment& transition_exp_used = targeted_exp_subsampled; + + // Change the feature finding parameters: + // - no RT score (since we don't know the correct retention time) + // - no RT window + // - no elution model score + // - no peak quality (use all peaks) + // - if best peptides should be used, use peak quality + MRMFeatureFinderScoring featureFinder; + Param feature_finder_param(default_ffparam); + feature_finder_param.setValue("Scores:use_rt_score", "false"); + feature_finder_param.setValue("Scores:use_elution_model_score", "false"); + feature_finder_param.setValue("rt_extraction_window", -1.0); + feature_finder_param.setValue("stop_report_after_feature", 1); + feature_finder_param.setValue("TransitionGroupPicker:PeakPickerChromatogram:signal_to_noise", 1.0); // set to 1.0 in all cases + feature_finder_param.setValue("TransitionGroupPicker:compute_peak_quality", "false"); // no peak quality -> take all peaks! + feature_finder_param.setValue("TransitionGroupPicker:compute_peak_quality", "true"); + feature_finder_param.setValue("TransitionGroupPicker:minimal_quality", irt_detection_param.getValue("InitialQualityCutoff")); + featureFinder.setParameters(feature_finder_param); + featureFinder.setStrictFlag(false); // Since we do not know if features are correct, we should not be strict + + FeatureMap featureFile; // for results + OpenMS::MRMFeatureFinderScoring::TransitionGroupMapType transition_group_map; // for results + std::vector empty_swath_maps; + TransformationDescription empty_trafo; // empty transformation + + // Prepare the data with the chromatograms + boost::shared_ptr xic_map(new PeakMap); + xic_map->setChromatograms(chromatograms); + OpenSwath::SpectrumAccessPtr chromatogram_ptr = OpenSwath::SpectrumAccessPtr(new OpenMS::SpectrumAccessOpenMS(xic_map)); + + featureFinder.pickExperiment(chromatogram_ptr, featureFile, transition_exp_used, empty_trafo, empty_swath_maps, transition_group_map); + + // 4. Find most likely correct feature for each compound and add it to the + // "pairs" vector by computing pairs of iRT and real RT. + // + // Note that the quality threshold will only be applied if + // estimateBestPeptides is true + std::vector > pairs; // store the RT pairs to write the output trafoXML + std::map best_features = OpenSwathHelper::simpleFindBestFeature(transition_group_map, + true, irt_detection_param.getValue("OverallQualityCutoff")); + OPENMS_LOG_DEBUG << "Extracted best features: " << best_features.size() << std::endl; + + // Create pairs vector and store peaks + std::map trgrmap_allpeaks; // store all peaks above cutoff + for (std::map::iterator it = best_features.begin(); it != best_features.end(); ++it) + { + pairs.emplace_back(it->second, PeptideRTMap[it->first]); // pair + if (transition_group_map.find(it->first) != transition_group_map.end()) + { + trgrmap_allpeaks[ it->first ] = &transition_group_map[ it->first]; + } + } - // Change the feature finding parameters: - // - no RT score (since we don't know the correct retention time) - // - no RT window - // - no elution model score - // - no peak quality (use all peaks) - // - if best peptides should be used, use peak quality - MRMFeatureFinderScoring featureFinder; - Param feature_finder_param(default_ffparam); - feature_finder_param.setValue("Scores:use_rt_score", "false"); - feature_finder_param.setValue("Scores:use_elution_model_score", "false"); - feature_finder_param.setValue("rt_extraction_window", -1.0); - feature_finder_param.setValue("stop_report_after_feature", 1); - feature_finder_param.setValue("TransitionGroupPicker:PeakPickerChromatogram:signal_to_noise", 1.0); // set to 1.0 in all cases - feature_finder_param.setValue("TransitionGroupPicker:compute_peak_quality", "false"); // no peak quality -> take all peaks! - feature_finder_param.setValue("TransitionGroupPicker:compute_peak_quality", "true"); - feature_finder_param.setValue("TransitionGroupPicker:minimal_quality", irt_detection_param.getValue("InitialQualityCutoff")); - featureFinder.setParameters(feature_finder_param); - featureFinder.setStrictFlag(false); // Since we do not know if features are correct, we should not be strict + // 5. Perform the outlier detection + std::vector > pairs_corrected; + String outlier_method = irt_detection_param.getValue("outlierMethod").toString(); + if (outlier_method == "iter_residual" || outlier_method == "iter_jackknife") + { + pairs_corrected = MRMRTNormalizer::removeOutliersIterative(pairs, min_rsq, min_coverage, + irt_detection_param.getValue("useIterativeChauvenet").toBool(), outlier_method); + } + else if (outlier_method == "ransac") + { + // First, estimate of the maximum deviation from RT that is tolerated: + // Because 120 min gradient can have around 4 min elution shift, we use + // a default value of 3 % of the gradient to find upper RT threshold (3.6 min). + double pcnt_rt_threshold = irt_detection_param.getValue("RANSACMaxPercentRTThreshold"); + double max_rt_threshold = (RTRange.second - RTRange.first) * pcnt_rt_threshold / 100.0; + + pairs_corrected = MRMRTNormalizer::removeOutliersRANSAC(pairs, min_rsq, min_coverage, + irt_detection_param.getValue("RANSACMaxIterations"), max_rt_threshold, + irt_detection_param.getValue("RANSACSamplingSize")); + } + else if (outlier_method == "none") + { + pairs_corrected = pairs; + } + else + { + throw Exception::IllegalArgument(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, + String("Illegal argument '") + outlier_method + + "' used for outlierMethod (valid: 'iter_residual', 'iter_jackknife', 'ransac', 'none')."); + } + OPENMS_LOG_DEBUG << "Performed outlier detection, left with features: " << pairs_corrected.size() << std::endl; - FeatureMap featureFile; // for results - OpenMS::MRMFeatureFinderScoring::TransitionGroupMapType transition_group_map; // for results - std::vector empty_swath_maps; - TransformationDescription empty_trafo; // empty transformation + // 7. Select the "correct" peaks for m/z (and IM) correction (e.g. remove those not + // part of the linear regression) + std::map trgrmap_final; // store all peaks above cutoff + for (const auto& it : trgrmap_allpeaks) + { + if (it.second->getFeatures().empty() ) {continue;} + const MRMFeature& feat = it.second->getBestFeature(); + + // Check if the current feature is in the list of pairs used for the + // linear RT regression (using other features may result in wrong + // calibration values). + // Matching only by RT is not perfect but should work for most cases. + for (Size pit = 0; pit < pairs_corrected.size(); pit++) + { + if (fabs(feat.getRT() - pairs_corrected[pit].first ) < 1e-2) + { + trgrmap_final[ it.first ] = it.second; + break; + } + } + } - // Prepare the data with the chromatograms - boost::shared_ptr xic_map(new PeakMap); - xic_map->setChromatograms(chromatograms); - OpenSwath::SpectrumAccessPtr chromatogram_ptr = OpenSwath::SpectrumAccessPtr(new OpenMS::SpectrumAccessOpenMS(xic_map)); + // 8. Correct m/z (and IM) deviations using SwathMapMassCorrection + // m/z correction is done with the -irt_im_extraction parameters + SwathMapMassCorrection mc; + mc.setParameters(calibration_param); + + mc.correctMZ(trgrmap_final, targeted_exp, swath_maps, pasef); + mc.correctIM(trgrmap_final, targeted_exp, swath_maps, pasef, im_trafo); + + // 9. store RT transformation, using the selected model + trafo_out.setDataPoints(pairs_corrected); + Param model_params; + model_params.setValue("symmetric_regression", "false"); + model_params.setValue("span", irt_detection_param.getValue("lowess:span")); + model_params.setValue("num_nodes", irt_detection_param.getValue("b_spline:num_nodes")); + String model_type = irt_detection_param.getValue("alignmentMethod").toString(); + trafo_out.fitModel(model_type, model_params); + + OPENMS_LOG_DEBUG << "Final RT mapping:" << std::endl; + for (Size i = 0; i < pairs_corrected.size(); i++) + { + OPENMS_LOG_DEBUG << pairs_corrected[i].first << " " << pairs_corrected[i].second << std::endl; + } - featureFinder.pickExperiment(chromatogram_ptr, featureFile, transition_exp_used, empty_trafo, empty_swath_maps, transition_group_map); + /////////////////////////////////////////////////////////////////////////// + // TODO incorperate this into transformation description function instead + ///////////////////////////////////////////////////////////// - // 4. Find most likely correct feature for each compound and add it to the - // "pairs" vector by computing pairs of iRT and real RT. - // - // Note that the quality threshold will only be applied if - // estimateBestPeptides is true - std::vector > pairs; // store the RT pairs to write the output trafoXML - std::map best_features = OpenSwathHelper::simpleFindBestFeature(transition_group_map, - true, irt_detection_param.getValue("OverallQualityCutoff")); - OPENMS_LOG_DEBUG << "Extracted best features: " << best_features.size() << std::endl; + std::vector predicted_values(pairs_corrected.size()); // these are the predicted values by applying the trafo on the values + std::vector delta_true_predicted(pairs_corrected.size()); // these are the differences between the predicted and the true values + double sum_rt_differences = 0.0; + TransformationDescription trafo_inverse = trafo_out; + trafo_inverse.invert(); + for (size_t i = 0; i < pairs_corrected.size(); ++i) + { + predicted_values[i] = trafo_out.apply(pairs_corrected[i].first); + delta_true_predicted[i] = pairs_corrected[i].first - predicted_values[i]; + sum_rt_differences += delta_true_predicted[i]; // sum for computing the mean + OPENMS_LOG_DEBUG << "True - Predicted value: " << pairs_corrected[i].second << " " << predicted_values[i] << std::endl; + } + double mean_rt_difference = sum_rt_differences / delta_true_predicted.size(); - // Create pairs vector and store peaks - std::map trgrmap_allpeaks; // store all peaks above cutoff - for (std::map::iterator it = best_features.begin(); it != best_features.end(); ++it) - { - pairs.emplace_back(it->second, PeptideRTMap[it->first]); // pair - if (transition_group_map.find(it->first) != transition_group_map.end()) + // compute the variance + double rt_variance = 0.0; + for (double delta: delta_true_predicted) { - trgrmap_allpeaks[ it->first ] = &transition_group_map[ it->first]; + rt_variance += (delta - mean_rt_difference) * (delta - mean_rt_difference); } - } + rt_variance /= delta_true_predicted.size(); + double rt_stdev = sqrt(rt_variance); + double rt_extraction_window = rt_stdev * 3.0 * 2; // 3 times the standard deviation, doubled since the window is the full length (not half) + std::cout << "RT variance: " << rt_stdev << std::endl; + std::cout << "RT Extraction Window to Use: " << rt_extraction_window << std::endl; + - // 5. Perform the outlier detection - std::vector > pairs_corrected; - String outlier_method = irt_detection_param.getValue("outlierMethod").toString(); - if (outlier_method == "iter_residual" || outlier_method == "iter_jackknife") - { - pairs_corrected = MRMRTNormalizer::removeOutliersIterative(pairs, min_rsq, min_coverage, - irt_detection_param.getValue("useIterativeChauvenet").toBool(), outlier_method); - } - else if (outlier_method == "ransac") - { - // First, estimate of the maximum deviation from RT that is tolerated: - // Because 120 min gradient can have around 4 min elution shift, we use - // a default value of 3 % of the gradient to find upper RT threshold (3.6 min). - double pcnt_rt_threshold = irt_detection_param.getValue("RANSACMaxPercentRTThreshold"); - double max_rt_threshold = (RTRange.second - RTRange.first) * pcnt_rt_threshold / 100.0; + std::cout << "End of doDataNormalization_ method" << std::endl; - pairs_corrected = MRMRTNormalizer::removeOutliersRANSAC(pairs, min_rsq, min_coverage, - irt_detection_param.getValue("RANSACMaxIterations"), max_rt_threshold, - irt_detection_param.getValue("RANSACSamplingSize")); - } - else if (outlier_method == "none") - { - pairs_corrected = pairs; + trafo_out.printSummary(std::cout); } - else { - throw Exception::IllegalArgument(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, - String("Illegal argument '") + outlier_method + - "' used for outlierMethod (valid: 'iter_residual', 'iter_jackknife', 'ransac', 'none')."); - } - OPENMS_LOG_DEBUG << "Performed outlier detection, left with features: " << pairs_corrected.size() << std::endl; - // 7. Select the "correct" peaks for m/z (and IM) correction (e.g. remove those not - // part of the linear regression) - std::map trgrmap_final; // store all peaks above cutoff - for (const auto& it : trgrmap_allpeaks) - { - if (it.second->getFeatures().empty() ) {continue;} - const MRMFeature& feat = it.second->getBestFeature(); + ////////////////////////////////////////////////////////////////////////////////////// + // PART 2: NONLINEAR CALIBRATION + ////////////////////////////////////////////////////////////////////////////////////// + int nrRTBins = (int) irt_detection_param.getValue("NrRTBins") * 10.0; // have 10X number of precursors as linear + int minPeptidesPerBin = (int) irt_detection_param.getValue("MinPeptidesPerBin") * 10; // have 10X number of precursors as linear + OpenSwath::LightTargetedExperiment targeted_exp_subsampled = OpenSwathHelper::subsampleLibrary(targeted_exp, + nrRTBins, + minPeptidesPerBin ); - // Check if the current feature is in the list of pairs used for the - // linear RT regression (using other features may result in wrong - // calibration values). - // Matching only by RT is not perfect but should work for most cases. - for (Size pit = 0; pit < pairs_corrected.size(); pit++) + + // 2. Store the peptide retention times in an intermediate map + std::map PeptideRTMap; + for (Size i = 0; i < targeted_exp_subsampled.getCompounds().size(); i++) { - if (fabs(feat.getRT() - pairs_corrected[pit].first ) < 1e-2) + PeptideRTMap[targeted_exp_subsampled.getCompounds()[i].id] = targeted_exp_subsampled.getCompounds()[i].rt; + } + + + // 2. Extract the chromatograms + std::vector< OpenMS::MSChromatogram > chromatograms; + TransformationDescription trafo; // dummy + cp_irt.rt_extraction_window = rt_extraction_window; + this->simpleExtractChromatograms_(swath_maps, targeted_exp_subsampled, chromatograms, trafo_out, cp_irt, pasef, load_into_memory); + + + OPENMS_LOG_DEBUG << "Extracted number of chromatograms from iRT files: " << chromatograms.size() << std::endl; + + // 2. For each peptide, extract and score chromatograms + const OpenSwath::LightTargetedExperiment& transition_exp_used = targeted_exp_subsampled; + + // Change the feature finding parameters: + // - no RT score (since we don't know the correct retention time) + // - no RT window + // - no elution model score + // - no peak quality (use all peaks) + // - if best peptides should be used, use peak quality + MRMFeatureFinderScoring featureFinder; + Param feature_finder_param(default_ffparam); + feature_finder_param.setValue("Scores:use_rt_score", "false"); + feature_finder_param.setValue("Scores:use_elution_model_score", "false"); + feature_finder_param.setValue("rt_extraction_window", -1.0); + feature_finder_param.setValue("stop_report_after_feature", 1); + feature_finder_param.setValue("TransitionGroupPicker:PeakPickerChromatogram:signal_to_noise", 1.0); // set to 1.0 in all cases + feature_finder_param.setValue("TransitionGroupPicker:compute_peak_quality", "false"); // no peak quality -> take all peaks! + feature_finder_param.setValue("TransitionGroupPicker:compute_peak_quality", "true"); + feature_finder_param.setValue("TransitionGroupPicker:minimal_quality", irt_detection_param.getValue("InitialQualityCutoff")); + featureFinder.setParameters(feature_finder_param); + featureFinder.setStrictFlag(false); // Since we do not know if features are correct, we should not be strict + + FeatureMap featureFile; // for results + OpenMS::MRMFeatureFinderScoring::TransitionGroupMapType transition_group_map; // for results + std::vector empty_swath_maps; + TransformationDescription empty_trafo; // empty transformation + + // Prepare the data with the chromatograms + boost::shared_ptr xic_map(new PeakMap); + xic_map->setChromatograms(chromatograms); + OpenSwath::SpectrumAccessPtr chromatogram_ptr = OpenSwath::SpectrumAccessPtr(new OpenMS::SpectrumAccessOpenMS(xic_map)); + + featureFinder.pickExperiment(chromatogram_ptr, featureFile, transition_exp_used, empty_trafo, empty_swath_maps, transition_group_map); + + // 4. Find most likely correct feature for each compound and add it to the + // "pairs" vector by computing pairs of iRT and real RT. + // + // Note that the quality threshold will only be applied if + // estimateBestPeptides is true + std::vector > pairs; // store the RT pairs to write the output trafoXML + std::map best_features = OpenSwathHelper::simpleFindBestFeature(transition_group_map, + true, irt_detection_param.getValue("OverallQualityCutoff")); + OPENMS_LOG_DEBUG << "Extracted best features: " << best_features.size() << std::endl; + + // Create pairs vector and store peaks + std::map trgrmap_allpeaks; // store all peaks above cutoff + for (std::map::iterator it = best_features.begin(); it != best_features.end(); ++it) + { + pairs.emplace_back(it->second, PeptideRTMap[it->first]); // pair + if (transition_group_map.find(it->first) != transition_group_map.end()) { - trgrmap_final[ it.first ] = it.second; - break; + trgrmap_allpeaks[ it->first ] = &transition_group_map[ it->first]; } } - } - // 8. Correct m/z (and IM) deviations using SwathMapMassCorrection - // m/z correction is done with the -irt_im_extraction parameters - SwathMapMassCorrection mc; - mc.setParameters(calibration_param); + // 5. Perform the outlier detection + std::vector > pairs_corrected; + String outlier_method = irt_detection_param.getValue("outlierMethod").toString(); + if (outlier_method == "iter_residual" || outlier_method == "iter_jackknife") + { + pairs_corrected = MRMRTNormalizer::removeOutliersIterative(pairs, min_rsq, min_coverage, + irt_detection_param.getValue("useIterativeChauvenet").toBool(), outlier_method); + } + else if (outlier_method == "ransac") + { + // First, estimate of the maximum deviation from RT that is tolerated: + // Because 120 min gradient can have around 4 min elution shift, we use + // a default value of 3 % of the gradient to find upper RT threshold (3.6 min). + double pcnt_rt_threshold = irt_detection_param.getValue("RANSACMaxPercentRTThreshold"); + double max_rt_threshold = (RTRange.second - RTRange.first) * pcnt_rt_threshold / 100.0; + + pairs_corrected = MRMRTNormalizer::removeOutliersRANSAC(pairs, min_rsq, min_coverage, + irt_detection_param.getValue("RANSACMaxIterations"), max_rt_threshold, + irt_detection_param.getValue("RANSACSamplingSize")); + } + else if (outlier_method == "none") + { + pairs_corrected = pairs; + } + else + { + throw Exception::IllegalArgument(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, + String("Illegal argument '") + outlier_method + + "' used for outlierMethod (valid: 'iter_residual', 'iter_jackknife', 'ransac', 'none')."); + } + OPENMS_LOG_DEBUG << "Performed outlier detection, left with features: " << pairs_corrected.size() << std::endl; - mc.correctMZ(trgrmap_final, targeted_exp, swath_maps, pasef); - mc.correctIM(trgrmap_final, targeted_exp, swath_maps, pasef, im_trafo); + // 7. Select the "correct" peaks for m/z (and IM) correction (e.g. remove those not + // part of the linear regression) + std::map trgrmap_final; // store all peaks above cutoff + for (const auto& it : trgrmap_allpeaks) + { + if (it.second->getFeatures().empty() ) {continue;} + const MRMFeature& feat = it.second->getBestFeature(); + + // Check if the current feature is in the list of pairs used for the + // linear RT regression (using other features may result in wrong + // calibration values). + // Matching only by RT is not perfect but should work for most cases. + for (Size pit = 0; pit < pairs_corrected.size(); pit++) + { + if (fabs(feat.getRT() - pairs_corrected[pit].first ) < 1e-2) + { + trgrmap_final[ it.first ] = it.second; + break; + } + } + } - // 9. store RT transformation, using the selected model - TransformationDescription trafo_out; - trafo_out.setDataPoints(pairs_corrected); - Param model_params; - model_params.setValue("symmetric_regression", "false"); - model_params.setValue("span", irt_detection_param.getValue("lowess:span")); - model_params.setValue("num_nodes", irt_detection_param.getValue("b_spline:num_nodes")); - String model_type = irt_detection_param.getValue("alignmentMethod").toString(); - trafo_out.fitModel(model_type, model_params); + // 8. Correct m/z (and IM) deviations using SwathMapMassCorrection + // m/z correction is done with the -irt_im_extraction parameters + SwathMapMassCorrection mc; + mc.setParameters(calibration_param); + + mc.correctMZ(trgrmap_final, targeted_exp, swath_maps, pasef); + mc.correctIM(trgrmap_final, targeted_exp, swath_maps, pasef, im_trafo); + + // 9. store RT transformation, using the selected model + TransformationDescription trafo_out; + trafo_out.setDataPoints(pairs_corrected); + Param model_params; + model_params.setValue("symmetric_regression", "false"); + model_params.setValue("span", irt_detection_param.getValue("lowess:span")); + model_params.setValue("num_nodes", irt_detection_param.getValue("b_spline:num_nodes")); + model_params.setValue("alignmentMethod", "lowess"); + String model_type = irt_detection_param.getValue("alignmentMethod").toString(); + trafo_out.fitModel(model_type, model_params); + + OPENMS_LOG_DEBUG << "Final RT mapping:" << std::endl; + for (Size i = 0; i < pairs_corrected.size(); i++) + { + OPENMS_LOG_DEBUG << pairs_corrected[i].first << " " << pairs_corrected[i].second << std::endl; + } - OPENMS_LOG_DEBUG << "Final RT mapping:" << std::endl; - for (Size i = 0; i < pairs_corrected.size(); i++) - { - OPENMS_LOG_DEBUG << pairs_corrected[i].first << " " << pairs_corrected[i].second << std::endl; - } + /////////////////////////////////////////////////////////////////////////// + // TODO incorperate this into transformation description function instead + ///////////////////////////////////////////////////////////// - /////////////////////////////////////////////////////////////////////////// - // TODO incorperate this into transformation description function instead - ///////////////////////////////////////////////////////////// + std::vector predicted_values(pairs_corrected.size()); // these are the predicted values by applying the trafo on the values + std::vector delta_true_predicted(pairs_corrected.size()); // these are the differences between the predicted and the true values + double sum_rt_differences = 0.0; + TransformationDescription trafo_inverse = trafo_out; + trafo_inverse.invert(); - std::vector predicted_values(pairs_corrected.size()); // these are the predicted values by applying the trafo on the values - std::vector delta_true_predicted(pairs_corrected.size()); // these are the differences between the predicted and the true values - double sum_rt_differences = 0.0; - for (size_t i = 0; i < pairs_corrected.size(); ++i) - { - predicted_values[i] = trafo_out.apply(pairs_corrected[i].first); - delta_true_predicted[i] = pairs_corrected[i].second - predicted_values[i]; - sum_rt_differences += delta_true_predicted[i]; // sum for computing the mean - OPENMS_LOG_DEBUG << "True - Predicted value: " << pairs_corrected[i].second << " " << predicted_values[i] << std::endl; - } - double mean_rt_difference = sum_rt_differences / delta_true_predicted.size(); + // note pairs correct first is exp RT, pairs corrected second is iRT + for (size_t i = 0; i < pairs_corrected.size(); ++i) + { + predicted_values[i] = trafo_inverse.apply(pairs_corrected[i].second); + delta_true_predicted[i] = pairs_corrected[i].first - predicted_values[i]; + sum_rt_differences += delta_true_predicted[i]; // sum for computing the mean + OPENMS_LOG_DEBUG << "True - Predicted value: " << pairs_corrected[i].second << " " << predicted_values[i] << std::endl; + } + double mean_rt_difference = sum_rt_differences / delta_true_predicted.size(); - // compute the variance - double rt_variance = 0.0; - for (double delta: delta_true_predicted) - { - rt_variance += (delta - mean_rt_difference) * (delta - mean_rt_difference); - } - rt_variance /= delta_true_predicted.size(); - double rt_stdev = sqrt(rt_variance); - double rt_extraction_window = rt_stdev * 3.0 * 2; // 3 times the standard deviation, doubled since the window is the full length (not half) - OPENMS_LOG_DEBUG << "RT variance: " << rt_stdev << std::endl; - OPENMS_LOG_DEBUG << "RT Extraction Window to Use: " << rt_extraction_window << std::endl; - + // compute the variance + double rt_variance = 0.0; + for (double delta: delta_true_predicted) + { + rt_variance += (delta - mean_rt_difference) * (delta - mean_rt_difference); + } + rt_variance /= delta_true_predicted.size(); + double rt_stdev = sqrt(rt_variance); + double rt_extraction_window = rt_stdev * 3.0 * 2; // 3 times the standard deviation, doubled since the window is the full length (not half) + std::cout << "RT variance: " << rt_stdev << std::endl; + std::cout << "RT Extraction Window to Use: " << rt_extraction_window << std::endl; + - OPENMS_LOG_DEBUG << "End of doDataNormalization_ method" << std::endl; + OPENMS_LOG_DEBUG << "End of doDataNormalization_ method" << std::endl; - trafo_out.printSummary(std::cout); + trafo_out.printSummary(std::cout); - this->endProgress(); - return trafo_out; + return trafo_out; + } } TransformationDescription OpenSwathCalibrationWorkflow::doDataNormalization_( diff --git a/src/topp/OpenSwathWorkflow.cpp b/src/topp/OpenSwathWorkflow.cpp index 1dc90818719..dd29213da71 100644 --- a/src/topp/OpenSwathWorkflow.cpp +++ b/src/topp/OpenSwathWorkflow.cpp @@ -925,6 +925,12 @@ class TOPPOpenSwathWorkflow min_rsq, min_coverage, feature_finder_param, cp_irt, linear_irt, no_calibration, debug_level, pasef, load_into_memory); + + if (!irt_trafo_out.empty()) + { + FileHandler().storeTransformations(irt_trafo_out, trafo_rtnorm, {FileTypes::TRANSFORMATIONXML}); + } + } if (nonlinear_irt_tr_file.empty()) {