diff --git a/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/DIAHelper.h b/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/DIAHelper.h index 31135e8f108..f3b2106be2f 100644 --- a/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/DIAHelper.h +++ b/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/DIAHelper.h @@ -10,61 +10,74 @@ #include #include - +#include #include - namespace OpenMS { + struct RangeMZ; + struct RangeMobility; class TheoreticalSpectrumGenerator; namespace DIAHelpers { - /** @brief Helper functions for the DIA scoring of OpenSWATH */ ///@{ /** - @brief Integrate intensity in a spectrum from start to end - - This function will integrate the intensity in a spectrum between mz_start - and mz_end, returning the total intensity and an intensity-weighted m/z - value. + @brief Integrate intensity in a spectrum from in @p mz_range (and @p im_range if defined) + returning the intensity-weighted m/z and im values as well as the total intensity. - @note If there is no signal, mz will be set to -1 and intensity to 0 + @note If there is no signal, @p mz and @p im will be set to -1 and intensity to 0 @return Returns true if a signal was found (and false if no signal was found) */ - OPENMS_DLLAPI bool integrateWindow(const OpenSwath::SpectrumPtr& spectrum, double mz_start, - double mz_end, double& mz, double& intensity, bool centroided = false); + OPENMS_DLLAPI bool integrateWindow(const OpenSwath::SpectrumPtr& spectrum, + double& mz, double& im, double& intensity, const RangeMZ& mz_range, const RangeMobility& im_range, bool centroided = false); /** - @brief Integrate intensities in a spectrum from start to end + @brief Integrate intensity in SpectrumSequence in range @p mz_range (and @p im_range if defined) + returning the intensity-weighted m/z and im values as well as the total intensity. + + @note If there is no signal, @p mz and @p im will be set to -1 and intensity to 0 + @return Returns true if a signal was found (and false if no signal was found) + */ + OPENMS_DLLAPI bool integrateWindow(const SpectrumSequence& spectrum, + double& mz, double& im, double& intensity, const RangeMZ& mz_range, const RangeMobility& im_range, bool centroided = false); + + /** + @brief Integrate intensities in a spectrum in range @p im_range (if defined) for multiple windows. + @param windows_center is a vector of the center location of the windows. + @param width is the width of the windows across mz + @param im_range is the range of the IM dimension (if defined) + @param remove_zero is a flag indicating whether to remove zero intensity windows + + Returns: + @param[out] integrated_windows_intensity is a vector of the integrated intensity for each window + @param[out] integrated_windows_mz is a vector of the integrated intensity-weighted m/z for each window + @param[out] integrated_windows_im is a vector of the integrated intensity-weighted im for each window */ OPENMS_DLLAPI void integrateWindows(const OpenSwath::SpectrumPtr& spectrum, //!< [in] Spectrum const std::vector& windows_center, //!< [in] center location double width, std::vector& integrated_windows_intensity, std::vector& integrated_windows_mz, + std::vector& integrated_windows_im, + const RangeMobility& im_range, bool remove_zero = false); - /** - @brief Integrate intensity in an ion mobility spectrum from start to end - - This function will integrate the intensity in a spectrum between mz_start - and mz_end, returning the total intensity and an intensity-weighted drift - time value. - @note If there is no signal, mz will be set to -1 and intensity to 0 + /** + @brief Integrate intensities of a SpectrumSequence in range @p im_range for multiple windows. */ - OPENMS_DLLAPI void integrateDriftSpectrum(const OpenSwath::SpectrumPtr& spectrum, - double mz_start, - double mz_end, - double & im, - double & intensity, - double drift_start, - double drift_end); - + OPENMS_DLLAPI void integrateWindows(const SpectrumSequence& spectrum, //!< [in] Spectrum + const std::vector& windows_center, //!< [in] center location + double width, + std::vector& integrated_windows_intensity, + std::vector& integrated_windows_mz, + std::vector& integrated_windows_im, + const RangeMobility& im_range, + bool remove_zero = false); /** @brief Adjust left/right window based on window and whether its ppm or not */ @@ -136,8 +149,27 @@ namespace OpenMS OPENMS_DLLAPI void extractFirst(const std::vector >& peaks, std::vector& mass); /// extract second from vector of pairs OPENMS_DLLAPI void extractSecond(const std::vector >& peaks, std::vector& mass); - + + /** @brief optionally convert a DIA extraction window from ppm to m/z + @param dia_extraction_window - how wide the extraction window is total (can be in m/z or ppm) + @param ppm - whether the extraction window is in ppm or not + @return the extraction window in m/z + */ + OPENMS_DLLAPI RangeMZ createMZRangePPM(double mz_ref, double dia_extraction_window, const bool ppm); + + + /** + @brief Helper function for integrating a spectrum. + */ + OPENMS_DLLAPI void integrateWindow_(const OpenSwath::SpectrumPtr& spectrum, + double & mz, + double & im, + double & intensity, + const RangeMZ & mz_range, + const RangeMobility & im_range, + bool centroided); + } + ///}@ - } -} +} //namespace OpenMS diff --git a/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/DIAPrescoring.h b/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/DIAPrescoring.h index 1add379b897..934cc769c96 100644 --- a/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/DIAPrescoring.h +++ b/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/DIAPrescoring.h @@ -13,6 +13,7 @@ #include #include +#include namespace OpenMS { @@ -53,8 +54,9 @@ namespace OpenMS and compute manhattan distance and dotprod score between spectrum intensities and simulated spectrum. */ - void score(OpenSwath::SpectrumPtr spec, + void score(const SpectrumSequence& spec, const std::vector& lt, + const RangeMobility& im_range, double& dotprod, double& manhattan) const; @@ -63,7 +65,7 @@ namespace OpenMS the SpectrumAccessPtr for all transitions groups in the LightTargetedExperiment. */ void operator()(const OpenSwath::SpectrumAccessPtr& swath_ptr, - OpenSwath::LightTargetedExperiment& transition_exp_used, + OpenSwath::LightTargetedExperiment& transition_exp_used, const RangeMobility& range_im, OpenSwath::IDataFrameWriter* ivw) const; }; diff --git a/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/DIAScoring.h b/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/DIAScoring.h index 01507fde644..794f4576a20 100644 --- a/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/DIAScoring.h +++ b/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/DIAScoring.h @@ -16,6 +16,8 @@ #include #include +#include + namespace OpenMS { class TheoreticalSpectrumGenerator; @@ -82,15 +84,17 @@ namespace OpenMS //@{ /// Isotope scores, see class description void dia_isotope_scores(const std::vector& transitions, - SpectrumPtrType spectrum, + SpectrumSequence& spectrum, OpenSwath::IMRMFeature* mrmfeature, + const RangeMobility& im_range, double& isotope_corr, double& isotope_overlap) const; /// Massdiff scores, see class description void dia_massdiff_score(const std::vector& transitions, - const SpectrumPtrType& spectrum, + const SpectrumSequence& spectrum, const std::vector& normalized_library_intensity, + const RangeMobility& im_range, double& ppm_score, double& ppm_score_weighted, std::vector& diff_ppm) const; @@ -103,23 +107,23 @@ namespace OpenMS @param ppm_score Resulting score @return False if no signal was found (and no sensible score calculated), true otherwise */ - bool dia_ms1_massdiff_score(double precursor_mz, SpectrumPtrType spectrum, + bool dia_ms1_massdiff_score(double precursor_mz, const SpectrumSequence& spectrum, const RangeMobility& im_range, double& ppm_score) const; /// Precursor isotope scores for precursors (peptides and metabolites) - void dia_ms1_isotope_scores_averagine(double precursor_mz, const SpectrumPtrType& spectrum, - double& isotope_corr, double& isotope_overlap, int charge_state) const; - void dia_ms1_isotope_scores(double precursor_mz, const SpectrumPtrType& spectrum, + void dia_ms1_isotope_scores_averagine(double precursor_mz, const SpectrumSequence& spectrum, int charge_state, RangeMobility& im_range, + double& isotope_corr, double& isotope_overlap) const; + void dia_ms1_isotope_scores(double precursor_mz, const std::vector& spectrum, RangeMobility& im_range, double& isotope_corr, double& isotope_overlap, const EmpiricalFormula& sum_formula) const; - /// b/y ion scores - void dia_by_ion_score(const SpectrumPtrType& spectrum, AASequence& sequence, - int charge, double& bseries_score, double& yseries_score) const; + void dia_by_ion_score(const SpectrumSequence& spectrum, AASequence& sequence, + int charge, const RangeMobility& im_range, double& bseries_score, double& yseries_score) const; /// Dotproduct / Manhattan score with theoretical spectrum - void score_with_isotopes(SpectrumPtrType spectrum, + void score_with_isotopes(SpectrumSequence& spectrum, const std::vector& transitions, + const RangeMobility& im_range, double& dotprod, double& manhattan) const; //@} @@ -137,8 +141,9 @@ namespace OpenMS /// Subfunction of dia_isotope_scores void diaIsotopeScoresSub_(const std::vector& transitions, - const SpectrumPtrType& spectrum, + const SpectrumSequence& spectrum, std::map& intensities, + const RangeMobility& im_range, double& isotope_corr, double& isotope_overlap) const; @@ -152,8 +157,8 @@ namespace OpenMS private: /** - @brief Determine whether the current m/z value is a monoisotopic peak - + @brief Determine whether the current m/z value is a monoisotopic peak + This function will try to determine whether the current peak is a monoisotopic peak or not. It will do so by searching for an intense peak at a lower m/z that could explain the current peak as part of a isotope @@ -166,7 +171,7 @@ namespace OpenMS @param nr_occurrences Will contain the maximum ratio of a peaks intensity compared to the monoisotopic peak intensity how often a peak is found at lower m/z than mono_mz with an intensity higher than mono_int. Multiple charge states are tested, see class parameter dia_nr_charges_ */ - void largePeaksBeforeFirstIsotope_(const SpectrumPtrType& spectrum, double mono_mz, double mono_int, int& nr_occurrences, double& max_ratio) const; + void largePeaksBeforeFirstIsotope_(const SpectrumSequence& spectrum, double mono_mz, double mono_int, int& nr_occurrences, double& max_ratio, const RangeMobility& im_range) const; /** @brief Compare an experimental isotope pattern to a theoretical one @@ -202,9 +207,8 @@ namespace OpenMS /// Get the intensities of isotopes around @p precursor_mz in experimental @p spectrum /// and fill @p isotopes_int. - void getIsotopeIntysFromExpSpec_(double precursor_mz, const SpectrumPtrType& spectrum, - std::vector& isotopes_int, - int charge_state) const; + void getIsotopeIntysFromExpSpec_(double precursor_mz, const SpectrumSequence& spectrum, int charge_state, const RangeMobility& im_range, + std::vector& isotopes_int) const; // Parameters double dia_extract_window_; @@ -219,4 +223,3 @@ namespace OpenMS TheoreticalSpectrumGenerator * generator; }; } - diff --git a/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/IonMobilityScoring.h b/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/IonMobilityScoring.h index 1095d98150c..3b714c810c8 100644 --- a/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/IonMobilityScoring.h +++ b/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/IonMobilityScoring.h @@ -14,17 +14,18 @@ #include #include #include - #include // scoring #include #include - namespace OpenMS { + class RangeMobility; + class RangeMZ; + /** @brief A class that calls the ion mobility scoring routines * * Use this class to invoke the individual OpenSWATH ion mobility scoring @@ -41,6 +42,18 @@ namespace OpenMS typedef OpenSwath::LightCompound CompoundType; typedef OpenSwath::LightTransition TransitionType; + struct MobilityPeak + { + double im; + double intensity; + MobilityPeak (); + MobilityPeak (double im_, double int_) : + im(im_), + intensity(int_) + {} + }; + typedef std::vector< MobilityPeak > IonMobilogram; + public: /// Constructor @@ -52,12 +65,11 @@ namespace OpenMS /** @brief Performs scoring of the ion mobility dimension in MS2 - @param spectrum The DIA MS2 spectrum found at the peak apex + @param spectrum sequence of segments of the DIA MS2 spectrum found at (and around) the peak apex @param transitions The transitions used for scoring @param scores The output scores - @param drift_lower Ion Mobility extraction start - @param drift_upper Ion Mobility extraction end @param drift_target Ion Mobility extraction target + @param im_range Ion Mobility extraction range @param dia_extraction_window_ m/z extraction width @param dia_extraction_ppm_ Whether m/z extraction width is in ppm @param use_spline Whether to use spline for fitting @@ -66,12 +78,11 @@ namespace OpenMS @return Populates additional scores in the @p scores object */ - static void driftScoring(const OpenSwath::SpectrumPtr& spectrum, + static void driftScoring(const SpectrumSequence& spectra, const std::vector & transitions, OpenSwath_Scores & scores, - const double drift_lower, - const double drift_upper, const double drift_target, + RangeMobility im_range, const double dia_extraction_window_, const bool dia_extraction_ppm_, const bool use_spline, @@ -80,11 +91,10 @@ namespace OpenMS /** @brief Performs scoring of the ion mobility dimension in MS1 - @param spectrum The DIA MS1 spectrum found at the peak apex + @param spectra vector containing the DIA MS1 spectra found at (or around) the peak apex @param transitions The transitions used for scoring @param scores The output scores - @param drift_lower Ion Mobility extraction start - @param drift_upper Ion Mobility extraction end + @param im_range Ion Mobility extraction range @param drift_target Ion Mobility extraction target @param dia_extraction_window_ m/z extraction width @param dia_extraction_ppm_ Whether m/z extraction width is in ppm @@ -94,12 +104,11 @@ namespace OpenMS @return Populates additional scores in the @p scores object */ - static void driftScoringMS1(const OpenSwath::SpectrumPtr& spectrum, + static void driftScoringMS1(const SpectrumSequence& spectra, const std::vector & transitions, OpenSwath_Scores & scores, - const double drift_lower, - const double drift_upper, const double drift_target, + RangeMobility im_range, const double dia_extract_window_, const bool dia_extraction_ppm_, const bool use_spline, @@ -108,13 +117,12 @@ namespace OpenMS /** @brief Performs scoring of the ion mobility dimension in MS1 and MS2 (contrast) - @param spectrum The DIA MS2 spectrum found at the peak apex - @param ms1spectrum The DIA MS1 spectrum found at the peak apex + @param spectra Vector of the DIA MS2 spectrum found in SpectrumSequence object (can contain 1 or multiple spectra centered around peak apex) + @param ms1spectrum The DIA MS1 spectrum found in SpectrumSequence object (can contain 1 or multiple spectra centered around peak apex) @param transitions The transitions used for scoring @param scores The output scores - @param drift_lower Ion Mobility extraction start - @param drift_upper Ion Mobility extraction end @param drift_target Ion Mobility extraction target + @param im_range the ion mobility range @param dia_extraction_window_ m/z extraction width @param dia_extraction_ppm_ Whether m/z extraction width is in ppm @param use_spline Whether to use spline for fitting @@ -123,14 +131,64 @@ namespace OpenMS @return Populates additional scores in the @p scores object */ - static void driftScoringMS1Contrast(const OpenSwath::SpectrumPtr& spectrum, const OpenSwath::SpectrumPtr& ms1spectrum, + static void driftScoringMS1Contrast(const SpectrumSequence& spectra, const SpectrumSequence& ms1spectrum, const std::vector & transitions, OpenSwath_Scores & scores, - const double drift_lower, - const double drift_upper, + RangeMobility im_range, const double dia_extract_window_, const bool dia_extraction_ppm_, const double drift_extra); + + /** + * @brief computes ion mobilogram to be used in scoring based on mz_range and im_range. + * Also integrates intensity in the resulting ion mobility mobilogram in mz_range and im_range across all the entire SpectrumSequence. + * @note If there is no signal, mz will be set to -1 and intensity to 0 + * @param[in] SpectrumSequence raw data in a spectrumSequence object (can contain 1 or multiple spectra centered around peak apex) + * @param[in] mz_range the range across mz to extract + * @param[in] im_range the range across im to extract + * @param[out] im computed weighted average ion mobility + * @param[out] intensity intensity computed intensity + * @param[out] res outputted ion mobilogram + * @param[in] eps minimum distance to allow for two seperate points + */ + static void computeIonMobilogram(const SpectrumSequence& spectra, + const RangeMZ & mz_range, + const RangeMobility & im_range, + double & im, + double & intensity, + IonMobilogram& res, + double eps); + + + private: + /** + * @brief helper function to computeIonMobilogram. Discretizes ion mobility values into a grid. + **/ + static std::vector computeGrid_(const std::vector< IonMobilogram >& mobilograms, double eps); + + + /* + @brief Extracts ion mobility values projected onto a grid + + For a given ion mobility profile and a grid, compute an ion mobilogram + across the grid for each ion mobility data point. Returns two data arrays + for the ion mobilogram: intensity (y) and ion mobility (x). Zero values are + inserted if no data point was found for a given grid value. + + @param profile The ion mobility data + @param im_grid The grid to be used + @param al_int_values The intensity vector (y) + @param al_im_values The ion mobility vector (x) + @param eps Epsilon used for computing the ion mobility grid + @param max_peak_idx The grid position of the maximum + */ + static void alignToGrid_(const IonMobilogram& profile, + const std::vector& im_grid, + std::vector< double >& al_int_values, + std::vector< double >& al_im_values, + double eps, + Size & max_peak_idx); + }; } diff --git a/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/MRMFeatureFinderScoring.h b/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/MRMFeatureFinderScoring.h index 1f3679408c1..c220fea6845 100644 --- a/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/MRMFeatureFinderScoring.h +++ b/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/MRMFeatureFinderScoring.h @@ -263,6 +263,7 @@ namespace OpenMS int stop_report_after_feature_; bool write_convex_hull_; bool strict_; + bool use_ms1_ion_mobility_; String scoring_model_; // scoring parameters diff --git a/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/OpenSwathScores.h b/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/OpenSwathScores.h index b300af5a537..f59ee8b3a8d 100644 --- a/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/OpenSwathScores.h +++ b/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/OpenSwathScores.h @@ -75,7 +75,7 @@ namespace OpenMS double weighted_coelution_score = 0; double weighted_xcorr_shape = 0; double weighted_massdev_score = 0; - + double ms1_xcorr_coelution_score = -1; double ms1_xcorr_coelution_contrast_score = 0; double ms1_xcorr_coelution_combined_score = 0; @@ -148,9 +148,9 @@ namespace OpenMS double calculate_lda_prescore(const OpenSwath_Scores& scores) const; /** @brief A scoring model for peak groups with a single transition - * + * * Manually derived scoring model for single transition peakgroups, only - * uses norm_rt_score, log_sn_score, and elution_model_fit_score. + * uses norm_rt_score, log_sn_score, and elution_model_fit_score. * * @returns A score which is better when more negative * diff --git a/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/OpenSwathScoring.h b/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/OpenSwathScoring.h index d7a89bdfad7..dc5b5e71c2b 100644 --- a/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/OpenSwathScoring.h +++ b/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/OpenSwathScoring.h @@ -23,24 +23,38 @@ #include #include +//logging +#include + +class RangeMZ; +class RangeMobility; + namespace OpenMS { /** @brief A class that calls the scoring routines * * Use this class to invoke the individual OpenSWATH scoring routines. - * + * */ - class OPENMS_DLLAPI OpenSwathScoring + class OPENMS_DLLAPI OpenSwathScoring { typedef OpenSwath::LightCompound CompoundType; typedef OpenSwath::LightTransition TransitionType; + enum class SpectrumAdditionMethod + { + ADDITION, + RESAMPLE + }; + double rt_normalization_factor_; double spacing_for_spectra_resampling_; int add_up_spectra_; - std::string spectra_addition_method_; + SpectrumAdditionMethod spectra_addition_method_; double im_drift_extra_pcnt_; OpenSwath_Scores_Usage su_; + bool use_ms1_ion_mobility_; // whether to use MS1 ion mobility extraction in DIA scores + const std::string ION_MOBILITY_DESCRIPTION = "Ion Mobility"; public: @@ -66,7 +80,8 @@ namespace OpenMS double spacing_for_spectra_resampling, const double drift_extra, const OpenSwath_Scores_Usage & su, - const std::string& spectrum_addition_method); + const std::string& spectrum_addition_method, + bool use_ms1_ion_mobility); /** @brief Score a single peakgroup in a chromatogram using only chromatographic properties. * @@ -93,7 +108,7 @@ namespace OpenMS std::vector& signal_noise_estimators, OpenSwath_Scores & scores) const; - /** @brief Score identification transitions against detection transitions of a single peakgroup + /** @brief Score identification transitions against detection transitions of a single peakgroup * in a chromatogram using only chromatographic properties. * * This function only uses the chromatographic properties (coelution, @@ -123,7 +138,7 @@ namespace OpenMS * peptide object. Both contain information about the expected elution time * on the chromatography and the relative intensity of the transitions. * - * The scores are returned in the OpenSwath_Scores object. + * The scores are returned in the OpenSwath_Scores object. * * @param imrmfeature The feature to be scored * @param transitions The library transition to score the feature against @@ -140,7 +155,7 @@ namespace OpenMS /** @brief Score a single chromatographic feature using DIA / SWATH scores. * - * The scores are returned in the OpenSwath_Scores object. + * The scores are returned in the OpenSwath_Scores object. * * @param imrmfeature The feature to be scored * @param transitions The library transition to score the feature against @@ -150,8 +165,8 @@ namespace OpenMS * @param pep The peptide corresponding to the library transitions * @param scores The object to store the result * @param mzerror_ppm m/z and mass error (in ppm) for all transitions - * @param drift_lower Drift time lower extraction boundary - * @param drift_upper Drift time upper extraction boundary + * @param[in] drift_target target drift value + * @param[in] range_im drift time lower and upper bounds * */ void calculateDIAScores(OpenSwath::IMRMFeature* imrmfeature, @@ -162,52 +177,48 @@ namespace OpenMS const CompoundType& compound, OpenSwath_Scores& scores, std::vector& mzerror_ppm, - const double drift_lower, - const double drift_upper, - const double drift_target); + const double drift_target, + const RangeMobility& range_im); /** @brief Score a single chromatographic feature using the precursor map. * - * The scores are returned in the OpenSwath_Scores object. + * The scores are returned in the OpenSwath_Scores object. * * @param ms1_map The MS1 (precursor ion map) from which the precursor spectra can be retrieved * @param diascoring DIA Scoring object to use for scoring * @param precursor_mz The m/z ratio of the precursor * @param rt The compound retention time + * @param compound the compound sequence + * @param im_range drift time lower and upper bounds * @param scores The object to store the result - * @param drift_lower Drift time lower extraction boundary - * @param drift_upper Drift time upper extraction boundary * */ - void calculatePrecursorDIAScores(const OpenSwath::SpectrumAccessPtr& ms1_map, + void calculatePrecursorDIAScores(const OpenSwath::SpectrumAccessPtr& ms1_map, const OpenMS::DIAScoring& diascoring, - double precursor_mz, - double rt, - const CompoundType& compound, - OpenSwath_Scores& scores, - double drift_lower, - double drift_upper); + double precursor_mz, + double rt, + const CompoundType& compound, + RangeMobility im_range, + OpenSwath_Scores& scores); /** @brief Score a single chromatographic feature using DIA / SWATH scores. * - * The scores are returned in the OpenSwath_Scores object. + * The scores are returned in the OpenSwath_Scores object. * * @param imrmfeature The feature to be scored * @param transitions The library transition to score the feature against * @param swath_maps The SWATH-MS (DIA) maps from which to retrieve full MS/MS spectra at the chromatographic peak apices + * @param range_im drift time lower and upper bounds * @param diascoring DIA Scoring object to use for scoring * @param scores The object to store the result - * @param drift_lower Drift time lower extraction boundary - * @param drift_upper Drift time upper extraction boundary * */ void calculateDIAIdScores(OpenSwath::IMRMFeature* imrmfeature, const TransitionType & transition, const std::vector& swath_maps, + RangeMobility& range_im, const OpenMS::DIAScoring & diascoring, - OpenSwath_Scores & scores, - double drift_lower, - double drift_upper); + OpenSwath_Scores & scores); /** @brief Computing the normalized library intensities from the transition objects * @@ -220,71 +231,57 @@ namespace OpenMS void getNormalized_library_intensities_(const std::vector & transitions, std::vector& normalized_library_intensity); - /** @brief Prepares a spectrum for DIA analysis (multiple map) - * - * This function will sum up (add) the intensities of multiple spectra from - * multiple swath maps (assuming these are SONAR maps of shifted precursor - * isolation windows) around the given retention time and return an - * "averaged" spectrum which may contain less noise. + /** @brief Prepares a spectrum for DIA analysis (single map) * - * @param[in] swath_maps The map(s) containing the spectra - * @param[in] RT The target retention time - * @param[in] nr_spectra_to_add How many spectra to add up - * @param drift_lower Drift time lower extraction boundary - * @param drift_upper Drift time upper extraction boundary + * This function will fetch a vector of spectrum pointers to be used in DIA analysis. + * If nr_spectra_to_add == 1, then a vector of length 1 will be returned * - * @return Added up spectrum + * Case #1: Non SONAR data and "simple" addition selected - Array of length "nr_spectra_to_add" returned corresponding with "nr_spectra_to_add" spectra + * Case #2: Non SONAR data and "resampling addition selected - Array of length 1 of the resampled spectrum returned + * Case #3: SONAR data - Array of length 1 containing the added/resampled spectrum returned * - */ - OpenSwath::SpectrumPtr fetchSpectrumSwath(std::vector swath_maps, - double RT, - int nr_spectra_to_add, - const double drift_lower, - const double drift_upper); - - /** @brief Prepares a spectrum for DIA analysis (single map) - * - * This function will sum up (add) the intensities of multiple spectra a single - * swath map (assuming these are regular SWATH / DIA maps) around the given + * For case #2 and #3 result is + * all spectra summed up (add) with the intensities of multiple spectra a single + * swath map (assuming these are regular SWATH / DIA maps) around the given * retention time and return an "averaged" spectrum which may contain less noise. * + * For case #1 this processing is done downstream in DIA scores to speed up computation time + * * @param[in] swath_map The map containing the spectra * @param[in] RT The target retention time * @param[in] nr_spectra_to_add How many spectra to add up - * @param drift_lower Drift time lower extraction boundary - * @param drift_upper Drift time upper extraction boundary - * - * @return Added up spectrum + * @param[in] range_im drift time lower and upper bounds + * @return Vector of spectra to be used * */ - OpenSwath::SpectrumPtr fetchSpectrumSwath(OpenSwath::SpectrumAccessPtr swath_map, - double RT, - int nr_spectra_to_add, - const double drift_lower, - const double drift_upper); + SpectrumSequence fetchSpectrumSwath(std::vector swath_maps, double RT, int nr_spectra_to_add, const RangeMobility& im_range); - protected: - /** @brief Returns an averaged spectrum + /** @brief Prepares a spectrum for DIA analysis (multiple map) + * + * This function will fetch a SpectrumSequence to be used in DIA analysis. + * If nr_spectra_to_add == 1, then a vector of length 1 will be returned. + * Spectra are prepared differently based on the condition + * Case #1: Non SONAR data and "simple" addition selected - Array of length "nr_spectra_to_add" returned corresponding with "nr_spectra_to_add" spectra + * Case #2: Non SONAR data and "resampling addition selected - Array of length 1 of the resampled spectrum returned + * Case #3: SONAR data - Array of length 1 containing the added/resampled spectrum returned + * + * For case #2 and #3 result is + * all spectra summed up (add) with the intensities of multiple spectra a single + * swath map (assuming these are regular SWATH / DIA maps) around the given + * retention time and return an "averaged" spectrum which may contain less noise. + * Spectra are also filtered and summed across drift time to transform an ion mobility spectrum into a non ion mobility spectrum * - * This function will sum up (add) the intensities of multiple spectra - * around the given retention time and return an "averaged" spectrum which - * may contain less noise. + * For case #1 this processing is done downstream in DIA scores to speed up computation time, furthermore drift time filtering is done downstream (these parameters are ignored) * * @param[in] swath_map The map containing the spectra * @param[in] RT The target retention time * @param[in] nr_spectra_to_add How many spectra to add up - * @param drift_lower Drift time lower extraction boundary - * @param drift_upper Drift time upper extraction boundary + * @param[in] im_range mobility range, only used if resampling spectrum addition chosen + * + * @return Vector of spectra to be used * - * @return Added up spectrum */ - OpenSwath::SpectrumPtr getAddedSpectra_(const OpenSwath::SpectrumAccessPtr& swath_map, - double RT, - int nr_spectra_to_add, - const double drift_lower, - const double drift_upper); - + SpectrumSequence fetchSpectrumSwath(OpenSwath::SpectrumAccessPtr swath_map, double RT, int nr_spectra_to_add, const RangeMobility& im_range); }; -} - +} \ No newline at end of file diff --git a/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/OpenSwathWorkflow.h b/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/OpenSwathWorkflow.h index 7d8089cdf50..dcd2447376e 100644 --- a/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/OpenSwathWorkflow.h +++ b/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/OpenSwathWorkflow.h @@ -215,7 +215,7 @@ namespace OpenMS * calibrants) perform RT and m/z correction in SWATH-MS data. Currently * supports (non-)linear correction of RT against library RT as well * as (non-)linear correction of m/z error as a function of m/z. - * + * * @note The relevant algorithms are implemented in MRMRTNormalizer for RT * calibration and SwathMapMassCorrection for m/z calibration. * @@ -356,10 +356,10 @@ namespace OpenMS /** * @brief Execute all steps in an \ref TOPP_OpenSwathWorkflow "OpenSwath" analysis * - * The workflow will perform a complete OpenSWATH analysis. Optionally, - * a calibration of m/z and retention time (mapping peptides to normalized - * space and correcting m/z error) can be performed beforehand using the - * OpenSwathCalibrationWorkflow class. + * The workflow will perform a complete OpenSWATH analysis. Optionally, + * a calibration of m/z and retention time (mapping peptides to normalized + * space and correcting m/z error) can be performed beforehand using the + * OpenSwathCalibrationWorkflow class. * * For diaPASEF workflows where ion mobility windows are overlapping, precursors may be found in multiple SWATHs. * In this case, precursors are only extracted from the SWATH in which they are most centered across ion mobility diff --git a/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/SpectrumAddition.h b/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/SpectrumAddition.h index 611a45f52f5..7b7bd0a5849 100644 --- a/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/SpectrumAddition.h +++ b/src/openms/include/OpenMS/ANALYSIS/OPENSWATH/SpectrumAddition.h @@ -16,8 +16,8 @@ namespace OpenMS { /** - @brief The SpectrumAddition is used to add up individual spectra - + @brief The SpectrumAddition is used to add up individual spectra + It uses the given sampling rate to resample the spectra in m/z domain and then add them up. This may lead to a certain inaccuracy, especially if a inappropriate resampling rate is chosen. @@ -29,15 +29,28 @@ namespace OpenMS public: /// adds up a list of Spectra by resampling them and then addition of intensities - static OpenSwath::SpectrumPtr addUpSpectra(const std::vector& all_spectra, + static OpenSwath::SpectrumPtr addUpSpectra(const SpectrumSequence& all_spectra, double sampling_rate, bool filter_zeros); + + /// adds up a list of ion mobility enhacned Spectra by resampling them and then addition of intensities. Currently this involves filtering to the desired IM extracion window and then performing addition across m/z and intensity. + static OpenSwath::SpectrumPtr addUpSpectra(const SpectrumSequence& all_spectra, + const RangeMobility& im_range, + double sampling_rate, + bool filter_zeros); + + /// Concatenates a spectrum sequence into a single spectrum. Values are sorted by m/z + static OpenSwath::SpectrumPtr concatenateSpectra(const SpectrumSequence& all_spectra); + + /// adds up a list of Spectra by resampling them and then addition of intensities - static OpenMS::MSSpectrum addUpSpectra(const std::vector< OpenMS::MSSpectrum>& all_spectra, + static OpenMS::MSSpectrum addUpSpectra(const std::vector& all_spectra, double sampling_rate, bool filter_zeros); + // sorts a spectrumPtr object by mz + static void sortSpectrumByMZ(OpenSwath::Spectrum&); }; } diff --git a/src/openms/include/OpenMS/INTERFACES/ISpectrumAccess.h b/src/openms/include/OpenMS/INTERFACES/ISpectrumAccess.h index cfdd7d32a67..64b74694f88 100644 --- a/src/openms/include/OpenMS/INTERFACES/ISpectrumAccess.h +++ b/src/openms/include/OpenMS/INTERFACES/ISpectrumAccess.h @@ -103,4 +103,3 @@ namespace Interfaces } //end namespace Interfaces } //end namespace OpenMS - diff --git a/src/openms/include/OpenMS/KERNEL/RangeManager.h b/src/openms/include/OpenMS/KERNEL/RangeManager.h index 6d302ce09d8..d2303cb091d 100644 --- a/src/openms/include/OpenMS/KERNEL/RangeManager.h +++ b/src/openms/include/OpenMS/KERNEL/RangeManager.h @@ -39,7 +39,13 @@ namespace OpenMS public: /// C'tor: initialize with empty range RangeBase() = default; - + + /// Cutom C'tor which sets the range to a singular point + RangeBase(const double single) : + min_(single), max_(single) + { + } + /// Custom C'tor to set min and max /// @throws Exception::InvalidRange if min > max RangeBase(const double min, const double max) : @@ -81,7 +87,7 @@ namespace OpenMS /// is the range empty (i.e. min > max)? bool isEmpty() const { - return min_ > max_; + return min_ > max_; } /// is @p value within [min, max]? @@ -97,11 +103,11 @@ namespace OpenMS } /** @name Accessors for min and max - + We use accessors, to keep range consistent (i.e. ensure that min <= max) */ - ///@{ - + ///@{ + /// sets the minimum (and the maximum, if uninitialized) void setMin(const double min) { @@ -137,7 +143,7 @@ namespace OpenMS min_ = std::min(min_, other.min_); max_ = std::max(max_, other.max_); } - + /// extend the range such that it includes the given @p value void extend(const double value) { @@ -205,16 +211,14 @@ namespace OpenMS shift(sandbox.max_ - max_); } } - } - /** @brief Scale the range of the dimension by a @p factor. A factor > 1 increases the range; factor < 1 decreases it. Let d = max - min; then min = min - d*(factor-1)/2, i.e. scale(1.5) extends the range by 25% on each side. - + Scaling an empty range will not have any effect. @param factor The multiplier to increase the range by @@ -274,7 +278,7 @@ namespace OpenMS double min_ = std::numeric_limits::max(); double max_ = std::numeric_limits::lowest(); }; - + OPENMS_DLLAPI std::ostream& operator<<(std::ostream& out, const RangeBase& b); struct OPENMS_DLLAPI RangeRT : public RangeBase { @@ -285,7 +289,7 @@ namespace OpenMS using RangeBase::RangeBase; // inherit C'tors from base /** @name Accessors for min and max - + We use accessors, to keep range consistent (i.e. ensure that min <= max) */ ///@{ @@ -333,7 +337,7 @@ namespace OpenMS return RangeBase::contains(inner_range); } }; - + OPENMS_DLLAPI std::ostream& operator<<(std::ostream& out, const RangeRT& range); struct OPENMS_DLLAPI RangeMZ : public RangeBase @@ -345,7 +349,7 @@ namespace OpenMS using RangeBase::RangeBase; // inherit C'tors from base /** @name Accessors for min and max - + We use accessors, to keep range consistent (i.e. ensure that min <= max) */ ///@{ @@ -403,7 +407,7 @@ namespace OpenMS using RangeBase::RangeBase; // inherit C'tors from base /** @name Accessors for min and max - + We use accessors, to keep range consistent (i.e. ensure that min <= max) */ ///@{ @@ -461,7 +465,7 @@ namespace OpenMS using RangeBase::RangeBase; // inherit C'tors from base /** @name Accessors for min and max - + We use accessors, to keep range consistent (i.e. ensure that min <= max) */ ///@{ @@ -511,7 +515,7 @@ namespace OpenMS }; OPENMS_DLLAPI std::ostream& operator<<(std::ostream& out, const RangeMobility& range); - + /// Enum listing state of dimensions (RangeBases) enum class HasRangeType { @@ -590,7 +594,7 @@ namespace OpenMS throw Exception::InvalidRange(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION , "No assignment took place (no dimensions in common!);"); } return *this; - } + } /// extend all dimensions which overlap with @p rhs to contain the range of @p rhs /// Dimensions which are not contained in @p rhs are left untouched. @@ -683,7 +687,7 @@ namespace OpenMS } } - + /// Clamp min/max of all overlapping dimensions to min/max of @p rhs /// Dimensions which are not contained in @p rhs or where rhs is empty are left untouched. /// @param rhs Range to clamp to @@ -747,7 +751,7 @@ namespace OpenMS assert((r_base != nullptr) && "No base class has this MSDim!"); return *r_base; } - + /// is any/some/all dimension in this range populated? HasRangeType hasRange() const { @@ -864,7 +868,7 @@ namespace OpenMS return out; } - /// use this class as a base class for containers, e.g. MSSpectrum. It forces them to implement 'updateRanges()' as a common interface + /// use this class as a base class for containers, e.g. MSSpectrum. It forces them to implement 'updateRanges()' as a common interface /// and provides a 'getRange()' member which saves casting to a range type manually template class RangeManagerContainer diff --git a/src/openms/source/ANALYSIS/OPENSWATH/DIAHelper.cpp b/src/openms/source/ANALYSIS/OPENSWATH/DIAHelper.cpp index c313bb64a3f..c6568e029fa 100644 --- a/src/openms/source/ANALYSIS/OPENSWATH/DIAHelper.cpp +++ b/src/openms/source/ANALYSIS/OPENSWATH/DIAHelper.cpp @@ -13,9 +13,12 @@ #include #include +#include #include +#include + namespace OpenMS::DIAHelpers { @@ -33,142 +36,295 @@ namespace OpenMS::DIAHelpers left -= mz_extract_window / 2.0; right += mz_extract_window / 2.0; } + // If the left value is now < 0, this is invalid correct it to be 0 + if (left < 0) + { + left = 0; + } } + // Helper for integrate window returns the sum of all intensities, sum of all ion mobilities and sum of all mz + // no expensive division calls + // assumes mz, im and intensity should already be initiated. + void integrateWindow_(const OpenSwath::SpectrumPtr& spectrum, + double & mz, + double & im, + double & intensity, + const RangeMZ & mz_range, + const RangeMobility & im_range, + bool centroided) + { + OPENMS_PRECONDITION(spectrum != nullptr, "precondition: Spectrum cannot be nullptr"); + OPENMS_PRECONDITION(spectrum->getMZArray() != nullptr, "precondition: Cannot integrate if no m/z is available."); + //OPENMS_PRECONDITION(!spectrum->getMZArray()->data.empty(), " precondition: Warning: Cannot integrate if spectrum is empty"); // This is not a failure should check for this afterwards + OPENMS_PRECONDITION(std::adjacent_find(spectrum->getMZArray()->data.begin(), + spectrum->getMZArray()->data.end(), std::greater()) == spectrum->getMZArray()->data.end(), + "Precondition violated: m/z vector needs to be sorted!" ); + OPENMS_PRECONDITION(spectrum->getMZArray()->data.size() == spectrum->getIntensityArray()->data.size(), "precondition: MZ and Intensity array need to have the same length."); + + // ion mobility specific preconditions + //OPENMS_PRECONDITION((im_range.isEmpty()) && (spectrum->getDriftTimeArray() != nullptr), "precondition: Cannot integrate with drift time if no drift time is available."); This is not a failure can handle this + OPENMS_PRECONDITION((spectrum->getDriftTimeArray() == nullptr) || (spectrum->getDriftTimeArray()->data.empty()) || (spectrum->getMZArray()->data.size() == spectrum->getDriftTimeArray()->data.size()), "precondition: MZ and Drift Time array need to have the same length."); + OPENMS_PRECONDITION(!centroided, throw Exception::NotImplemented(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION)); + + if ( spectrum->getMZArray()->data.empty() ) + { + OPENMS_LOG_WARN << "Warning: Cannot integrate if spectrum is empty" << std::endl; + return; + } + + // if im_range is set, than integrate across dirft time + if (!im_range.isEmpty()) // if imRange supplied, integrate across IM + { + if (spectrum->getDriftTimeArray() == nullptr) + { + //throw Exception::MissingInformation(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "Cannot integrate with drift time if no drift time is available"); + throw Exception::MissingInformation(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "Cannot integrate with drift time if no drift time is available"); + } + } + + if (!centroided) + { + // get the weighted average for noncentroided data. + // TODO this is not optimal if there are two peaks in this window (e.g. if the window is too large) + auto mz_arr_end = spectrum->getMZArray()->data.end(); + auto int_it = spectrum->getIntensityArray()->data.begin(); + + // this assumes that the spectra are sorted! + auto mz_it = std::lower_bound(spectrum->getMZArray()->data.begin(), mz_arr_end, mz_range.getMin()); + + // also advance intensity iterator now + auto iterator_pos = std::distance(spectrum->getMZArray()->data.begin(), mz_it); + std::advance(int_it, iterator_pos); + + double mz_end = mz_range.getMax(); // store the maximum mz value in a double to minimize function calls + + if ( !im_range.isEmpty() ) // integrate across im as well + { + auto im_it = spectrum->getDriftTimeArray()->data.begin(); + + // also advance ion mobility iterator now + std::advance(im_it, iterator_pos); + + // Start iteration from mz start, end iteration when mz value is larger than mz_end, only store only storing ion mobility values that are in the range + //while ( (mz_it != mz_arr_end) && (*mz_it < mz_end) ) + while ( (mz_it != mz_arr_end) && (*mz_it < mz_end) ) + { + if (im_range.contains(*im_it)) + { + intensity += (*int_it); + im += (*int_it) * (*im_it); + mz += (*int_it) * (*mz_it); + } + ++mz_it; + ++int_it; + ++im_it; + } + } + else // where do not have IM + { + while ( mz_it != mz_arr_end && *mz_it < mz_end ) + { + intensity += (*int_it); + mz += (*int_it) * (*mz_it); + + ++mz_it; + ++int_it; + } + } + } + else + { + throw Exception::NotImplemented(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION); + } + } + + void integrateWindows(const OpenSwath::SpectrumPtr& spectrum, const std::vector & windowsCenter, double width, std::vector & integratedWindowsIntensity, std::vector & integratedWindowsMZ, + std::vector & integratedWindowsIm, + const RangeMobility & range_im, bool remZero) { std::vector::const_iterator beg = windowsCenter.begin(); std::vector::const_iterator end = windowsCenter.end(); - double mz, intensity; + double mz, intensity, im; for (; beg != end; ++beg) { - double left = *beg - width / 2.0; - double right = *beg + width / 2.0; - if (integrateWindow(spectrum, left, right, mz, intensity, false)) + // assemble RangeMZ object based on window + RangeMZ range_mz(*beg); + range_mz.minSpanIfSingular(width); + + if (integrateWindow(spectrum, mz, im, intensity, range_mz, range_im, false)) { integratedWindowsIntensity.push_back(intensity); integratedWindowsMZ.push_back(mz); + integratedWindowsIm.push_back(im); } else if (!remZero) { integratedWindowsIntensity.push_back(0.); integratedWindowsMZ.push_back(*beg); + if ( !range_im.isEmpty() ) + { + integratedWindowsIm.push_back( range_im.center() ); // average drift time + } + else + { + integratedWindowsIm.push_back(-1); + } } } } - void integrateDriftSpectrum(const OpenSwath::SpectrumPtr& spectrum, - double mz_start, - double mz_end, - double & im, - double & intensity, - double drift_start, - double drift_end) - { - OPENMS_PRECONDITION(spectrum->getDriftTimeArray() != nullptr, "Cannot filter by drift time if no drift time is available."); - OPENMS_PRECONDITION(spectrum->getMZArray()->data.size() == spectrum->getIntensityArray()->data.size(), "MZ and Intensity array need to have the same length."); - OPENMS_PRECONDITION(spectrum->getMZArray()->data.size() == spectrum->getDriftTimeArray()->data.size(), "MZ and Drift Time array need to have the same length."); - OPENMS_PRECONDITION(std::adjacent_find(spectrum->getMZArray()->data.begin(), - spectrum->getMZArray()->data.end(), std::greater()) == spectrum->getMZArray()->data.end(), - "Precondition violated: m/z vector needs to be sorted!" ) - - im = 0; - intensity = 0; - // get the weighted average for noncentroided data. - // TODO this is not optimal if there are two peaks in this window (e.g. if the window is too large) - auto mz_arr_end = spectrum->getMZArray()->data.end(); - auto int_it = spectrum->getIntensityArray()->data.begin(); - auto im_it = spectrum->getDriftTimeArray()->data.begin(); + void integrateWindows(const SpectrumSequence& spectra, + const std::vector & windowsCenter, + double width, + std::vector & integratedWindowsIntensity, + std::vector & integratedWindowsMZ, + std::vector & integratedWindowsIm, + const RangeMobility& range_im, + bool remZero) + { - // this assumes that the spectra are sorted! - auto mz_it = std::lower_bound(spectrum->getMZArray()->data.begin(), mz_arr_end, mz_start); - auto mz_it_end = std::lower_bound(mz_it, mz_arr_end, mz_end); + double mz(-1), intensity(0), im(-1); + if (windowsCenter.empty()) + { + throw Exception::MissingInformation(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "No windows supplied!"); + return; + } - // also advance intensity and ion mobility iterator now - auto iterator_pos = std::distance(spectrum->getMZArray()->data.begin(), mz_it); - std::advance(int_it, iterator_pos); - std::advance(im_it, iterator_pos); + if (spectra.empty()) + { + OPENMS_LOG_WARN << "Warning: no spectra provided" << std::endl; + return; + } - // Iterate from mz start to end, only storing ion mobility values that are in the range - for (; mz_it != mz_it_end; ++mz_it, ++int_it, ++im_it) + std::vector::const_iterator beg = windowsCenter.begin(); + std::vector::const_iterator end = windowsCenter.end(); + for (; beg != end; ++beg) { - if ( *im_it >= drift_start && *im_it <= drift_end) + //assemble rangeMZ object based on windows + RangeMZ range_mz(*beg); + range_mz.minSpanIfSingular(width); + + if (integrateWindow(spectra, mz, im, intensity, range_mz, range_im, false)) { - intensity += (*int_it); - im += (*int_it) * (*im_it); + integratedWindowsIntensity.push_back(intensity); + integratedWindowsMZ.push_back(mz); + integratedWindowsIm.push_back(im); + } + else if (!remZero) + { + integratedWindowsIntensity.push_back(0.); + integratedWindowsMZ.push_back(*beg); // push back center of window + if ( !range_im.isEmpty() ) + { + integratedWindowsIm.push_back(range_im.center()); // push back average drift + } + else + { + integratedWindowsIm.push_back(-1); + } } } + } + + bool integrateWindow(const OpenSwath::SpectrumPtr& spectrum, + double & mz, + double & im, + double & intensity, + const RangeMZ & range_mz, + const RangeMobility & range_im, + bool centroided) + { + + // initiate the values + mz = 0; + im = 0; + intensity = 0; + + integrateWindow_(spectrum, mz, im, intensity, range_mz, range_im, centroided); + // Post processing get the weighted average mz and im by dividing my intensity if (intensity > 0.) { - im /= intensity; + mz /= intensity; + + if ( !range_im.isEmpty() ) + { + im /= intensity; + } + else + { + im = -1; + } + return true; } else { im = -1; + mz = -1; intensity = 0; + return false; } - } - bool integrateWindow(const OpenSwath::SpectrumPtr& spectrum, - double mz_start, - double mz_end, - double & mz, - double & intensity, - bool centroided) + bool integrateWindow(const SpectrumSequence& spectra, + double & mz, + double & im, + double & intensity, + const RangeMZ & range_mz, + const RangeMobility & range_im, + bool centroided) { - OPENMS_PRECONDITION(spectrum->getMZArray()->data.size() == spectrum->getIntensityArray()->data.size(), "MZ and Intensity array need to have the same length."); - OPENMS_PRECONDITION(std::adjacent_find(spectrum->getMZArray()->data.begin(), - spectrum->getMZArray()->data.end(), std::greater()) == spectrum->getMZArray()->data.end(), - "Precondition violated: m/z vector needs to be sorted!" ) - + // initiate the values mz = 0; + im = 0; intensity = 0; - if (!centroided) - { - // get the weighted average for noncentroided data. - // TODO this is not optimal if there are two peaks in this window (e.g. if the window is too large) - auto mz_arr_end = spectrum->getMZArray()->data.end(); - auto int_it = spectrum->getIntensityArray()->data.begin(); - - // this assumes that the spectra are sorted! - auto mz_it = std::lower_bound(spectrum->getMZArray()->data.begin(), mz_arr_end, mz_start); - auto mz_it_end = std::lower_bound(mz_it, mz_arr_end, mz_end); - // also advance intensity iterator now - auto iterator_pos = std::distance(spectrum->getMZArray()->data.begin(), mz_it); - std::advance(int_it, iterator_pos); - - for (; mz_it != mz_it_end; ++mz_it, ++int_it) + if (!spectra.empty()) + { + for (const auto& s : spectra) { - intensity += (*int_it); - mz += (*int_it) * (*mz_it); + integrateWindow_(s, mz, im, intensity, range_mz, range_im, centroided); } + // Post processing get the weighted average mz and im by dividing my intensity if (intensity > 0.) { mz /= intensity; + if ( !range_im.isEmpty() ) + { + im /= intensity; + } + else // if no IM set to -1 + { + im = -1; + } return true; } else { + // if (intensity <= 0) + im = -1; mz = -1; intensity = 0; return false; } - } else { - // not implemented - throw Exception::NotImplemented(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION); + // if (all_spectra.empty()) + OPENMS_LOG_WARN << "Warning: no spectra provided" << std::endl; + im = -1; + mz = -1; + intensity = 0; + return false; } } @@ -355,6 +511,21 @@ namespace OpenMS::DIAHelpers } } + RangeMZ createMZRangePPM(const double mzRef, const double dia_extraction_window, const bool is_ppm) + { + RangeMZ rangeMZ(mzRef); + if (is_ppm) + { + double ppm = Math::ppmToMass(dia_extraction_window, mzRef); + rangeMZ.minSpanIfSingular(ppm); + } + else + { + rangeMZ.minSpanIfSingular(dia_extraction_window); + } + return rangeMZ; + } + //modify masses by charge void modifyMassesByCharge( const std::vector >& isotopeSpec, diff --git a/src/openms/source/ANALYSIS/OPENSWATH/DIAPrescoring.cpp b/src/openms/source/ANALYSIS/OPENSWATH/DIAPrescoring.cpp index b14a3efc124..465198615dd 100644 --- a/src/openms/source/ANALYSIS/OPENSWATH/DIAPrescoring.cpp +++ b/src/openms/source/ANALYSIS/OPENSWATH/DIAPrescoring.cpp @@ -8,7 +8,7 @@ #include -#include +//#include #include #include #include @@ -53,7 +53,7 @@ namespace OpenMS } void DiaPrescore::operator()(const OpenSwath::SpectrumAccessPtr& swath_ptr, - OpenSwath::LightTargetedExperiment& transition_exp_used, + OpenSwath::LightTargetedExperiment& transition_exp_used, const RangeMobility& im_range, OpenSwath::IDataFrameWriter* ivw) const { //getParams(); @@ -76,8 +76,9 @@ namespace OpenMS for (UInt i = 0; i < swath_ptr->getNrSpectra(); ++i) { - - OpenSwath::SpectrumPtr spec = swath_ptr->getSpectrumById(i); + OpenSwath::SpectrumPtr s = swath_ptr->getSpectrumById(i); + SpectrumSequence spec; + spec.push_back(s); OpenSwath::SpectrumMeta specmeta = swath_ptr->getSpectrumMetaById(i); std::cout << "Processing Spectrum " << i << "RT " << specmeta.RT << std::endl; @@ -94,7 +95,7 @@ namespace OpenMS double score1; double score2; //OpenSwath::LightPeptide pep; - score(spec, beg->second, score1, score2); + score(spec, beg->second, im_range, score1, score2); score1v.push_back(score1); score2v.push_back(score2); @@ -107,8 +108,9 @@ namespace OpenMS } //end of for loop over spectra } - void DiaPrescore::score(OpenSwath::SpectrumPtr spec, + void DiaPrescore::score(const SpectrumSequence& spec, const std::vector& lt, + const RangeMobility& im_range, double& dotprod, double& manhattan) const { @@ -157,8 +159,8 @@ namespace OpenMS std::vector mzTheor, intTheor; DIAHelpers::extractFirst(spectrumWIso, mzTheor); DIAHelpers::extractSecond(spectrumWIso, intTheor); - std::vector intExp, mzExp; - DIAHelpers::integrateWindows(std::move(spec), mzTheor, dia_extract_window_, intExp, mzExp); + std::vector intExp, mzExp, imExp; + DIAHelpers::integrateWindows(spec, mzTheor, dia_extract_window_, intExp, mzExp, imExp, im_range); std::transform(intExp.begin(), intExp.end(), intExp.begin(), [](double val){return std::sqrt(val);}); std::transform(intTheor.begin(), intTheor.end(), intTheor.begin(), [](double val){return std::sqrt(val);}); diff --git a/src/openms/source/ANALYSIS/OPENSWATH/DIAScoring.cpp b/src/openms/source/ANALYSIS/OPENSWATH/DIAScoring.cpp index 8cb09535916..0b594d4e2c4 100644 --- a/src/openms/source/ANALYSIS/OPENSWATH/DIAScoring.cpp +++ b/src/openms/source/ANALYSIS/OPENSWATH/DIAScoring.cpp @@ -77,7 +77,7 @@ namespace OpenMS // generator->setParameters(p); } - DIAScoring::~DIAScoring() + DIAScoring::~DIAScoring() { delete generator; } @@ -98,36 +98,36 @@ namespace OpenMS /////////////////////////////////////////////////////////////////////////// // DIA / SWATH scoring - void DIAScoring::dia_isotope_scores(const std::vector& transitions, SpectrumPtrType spectrum, - OpenSwath::IMRMFeature* mrmfeature, double& isotope_corr, double& isotope_overlap) const + void DIAScoring::dia_isotope_scores(const std::vector& transitions, std::vector& spectrum, + OpenSwath::IMRMFeature* mrmfeature, const RangeMobility& im_range, double& isotope_corr, double& isotope_overlap) const { isotope_corr = 0; isotope_overlap = 0; // first compute a map of relative intensities from the feature, then compute the score std::map intensities; getFirstIsotopeRelativeIntensities_(transitions, mrmfeature, intensities); - diaIsotopeScoresSub_(transitions, std::move(spectrum), intensities, isotope_corr, isotope_overlap); + diaIsotopeScoresSub_(transitions, spectrum, intensities, im_range, isotope_corr, isotope_overlap); } void DIAScoring::dia_massdiff_score(const std::vector& transitions, - const SpectrumPtrType& spectrum, + const SpectrumSequence& spectrum, const std::vector& normalized_library_intensity, + const RangeMobility& im_range, double& ppm_score, - double& ppm_score_weighted, + double& ppm_score_weighted, std::vector& diff_ppm) const { + // Calculate the difference of the theoretical mass and the actually measured mass ppm_score = 0; ppm_score_weighted = 0; diff_ppm.clear(); for (std::size_t k = 0; k < transitions.size(); k++) { const TransitionType& transition = transitions[k]; - // Calculate the difference of the theoretical mass and the actually measured mass - double left(transition.getProductMZ()), right(transition.getProductMZ()); - DIAHelpers::adjustExtractionWindow(right, left, dia_extract_window_, dia_extraction_ppm_); - double mz, intensity; - bool signalFound = DIAHelpers::integrateWindow(spectrum, left, right, mz, intensity, dia_centroided_); + RangeMZ mz_range = DIAHelpers::createMZRangePPM(transition.getProductMZ(), dia_extract_window_, dia_extraction_ppm_); + double mz, intensity, im; + bool signalFound = DIAHelpers::integrateWindow(spectrum, mz, im, intensity, mz_range, im_range, dia_centroided_); // Continue if no signal was found - we therefore don't make a statement // about the mass difference if no signal is present. if (!signalFound) @@ -149,22 +149,21 @@ namespace OpenMS ppm_score /= transitions.size(); } - bool DIAScoring::dia_ms1_massdiff_score(double precursor_mz, SpectrumPtrType spectrum, - double& ppm_score) const + bool DIAScoring::dia_ms1_massdiff_score(double precursor_mz, const SpectrumSequence& spectrum, + const RangeMobility& im_range, double& ppm_score) const { ppm_score = -1; - double mz, intensity; + double mz, intensity, im; { // Calculate the difference of the theoretical mass and the actually measured mass - double left(precursor_mz), right(precursor_mz); - DIAHelpers::adjustExtractionWindow(right, left, dia_extract_window_, dia_extraction_ppm_); - bool signalFound = DIAHelpers::integrateWindow(std::move(spectrum), left, right, mz, intensity, dia_centroided_); + RangeMZ mz_range = DIAHelpers::createMZRangePPM(precursor_mz, dia_extract_window_, dia_extraction_ppm_); + bool signalFound = DIAHelpers::integrateWindow(spectrum, mz, im, intensity, mz_range, im_range, dia_centroided_); // Catch if no signal was found and replace it with the most extreme // value. Otherwise, calculate the difference in ppm. if (!signalFound) { - ppm_score = (right - left) / precursor_mz * 1000000; + ppm_score = Math::getPPMAbs(precursor_mz + mz_range.getSpan(), precursor_mz); return false; } else @@ -176,14 +175,14 @@ namespace OpenMS } /// Precursor isotope scores - void DIAScoring::dia_ms1_isotope_scores(double precursor_mz, const SpectrumPtrType& spectrum, - double& isotope_corr, double& isotope_overlap, const EmpiricalFormula& sum_formula) const + void DIAScoring::dia_ms1_isotope_scores(double precursor_mz, const std::vector& spectrum, + RangeMobility& im_range, double& isotope_corr, double& isotope_overlap, const EmpiricalFormula& sum_formula) const { // although precursor_mz can be received from the empirical formula (if non-empty), the actual precursor could be // slightly different. And also for compounds, usually the neutral sum_formula without adducts is given. // Therefore calculate the isotopes based on the formula but place them at precursor_mz std::vector isotopes_int; - getIsotopeIntysFromExpSpec_(precursor_mz, spectrum, isotopes_int, sum_formula.getCharge()); + getIsotopeIntysFromExpSpec_(precursor_mz, spectrum, sum_formula.getCharge(), im_range, isotopes_int); double max_ratio = 0; int nr_occurrences = 0; @@ -191,32 +190,29 @@ namespace OpenMS // calculate the scores: // isotope correlation (forward) and the isotope overlap (backward) scores isotope_corr = scoreIsotopePattern_(isotopes_int, sum_formula); - largePeaksBeforeFirstIsotope_(spectrum, precursor_mz, isotopes_int[0], nr_occurrences, max_ratio); + largePeaksBeforeFirstIsotope_(spectrum, precursor_mz, isotopes_int[0], nr_occurrences, max_ratio, im_range); isotope_overlap = max_ratio; } - void DIAScoring::getIsotopeIntysFromExpSpec_(double precursor_mz, const SpectrumPtrType& spectrum, - std::vector& isotopes_int, - int charge_state) const + void DIAScoring::getIsotopeIntysFromExpSpec_(double precursor_mz, const SpectrumSequence& spectrum, int charge_state, const RangeMobility& im_range, + std::vector& isotopes_int) const { double abs_charge = std::fabs(static_cast(charge_state)); for (int iso = 0; iso <= dia_nr_isotopes_; ++iso) { - double left = precursor_mz + iso * C13C12_MASSDIFF_U / abs_charge; - double right = left; - DIAHelpers::adjustExtractionWindow(right, left, dia_extract_window_, dia_extraction_ppm_); - double mz, intensity; - DIAHelpers::integrateWindow(spectrum, left, right, mz, intensity, dia_centroided_); + RangeMZ mz_range = DIAHelpers::createMZRangePPM(precursor_mz + iso * C13C12_MASSDIFF_U / abs_charge, dia_extract_window_, dia_extraction_ppm_); + double mz, intensity, im; + + DIAHelpers::integrateWindow(spectrum, mz, im, intensity, mz_range, im_range, dia_centroided_); isotopes_int.push_back(intensity); } } - void DIAScoring::dia_ms1_isotope_scores_averagine(double precursor_mz, const SpectrumPtrType& spectrum, - double& isotope_corr, double& isotope_overlap, - int charge_state) const + void DIAScoring::dia_ms1_isotope_scores_averagine(double precursor_mz, const SpectrumSequence& spectrum, int charge_state, RangeMobility& im_range, + double& isotope_corr, double& isotope_overlap) const { std::vector exp_isotopes_int; - getIsotopeIntysFromExpSpec_(precursor_mz, spectrum, exp_isotopes_int, charge_state); + getIsotopeIntysFromExpSpec_(precursor_mz, spectrum, charge_state, im_range, exp_isotopes_int); CoarseIsotopePatternGenerator solver(dia_nr_isotopes_ + 1); // NOTE: this is a rough estimate of the neutral mz value since we would not know the charge carrier for negative ions IsotopeDistribution isotope_dist = solver.estimateFromPeptideWeight(std::fabs(precursor_mz * charge_state)); @@ -226,28 +222,26 @@ namespace OpenMS // calculate the scores: // isotope correlation (forward) and the isotope overlap (backward) scores isotope_corr = scoreIsotopePattern_(exp_isotopes_int, isotope_dist); - largePeaksBeforeFirstIsotope_(spectrum, precursor_mz, exp_isotopes_int[0], nr_occurrences, max_ratio); + largePeaksBeforeFirstIsotope_(spectrum, precursor_mz, exp_isotopes_int[0], nr_occurrences, max_ratio, im_range); isotope_overlap = max_ratio; } - void DIAScoring::dia_by_ion_score(const SpectrumPtrType& spectrum, - AASequence& sequence, int charge, double& bseries_score, + void DIAScoring::dia_by_ion_score(const SpectrumSequence& spectrum, + AASequence& sequence, int charge, const RangeMobility& im_range, double& bseries_score, double& yseries_score) const { bseries_score = 0; yseries_score = 0; OPENMS_PRECONDITION(charge > 0, "Charge is a positive integer"); // for peptides, charge should be positive - double mz, intensity, left, right; + double mz, intensity, im; std::vector yseries, bseries; OpenMS::DIAHelpers::getBYSeries(sequence, bseries, yseries, generator, charge); for (const auto& b_ion_mz : bseries) { - left = b_ion_mz; - right = b_ion_mz; - DIAHelpers::adjustExtractionWindow(right, left, dia_extract_window_, dia_extraction_ppm_); + RangeMZ mz_range = DIAHelpers::createMZRangePPM(b_ion_mz, dia_extract_window_, dia_extraction_ppm_); - bool signalFound = DIAHelpers::integrateWindow(spectrum, left, right, mz, intensity, dia_centroided_); + bool signalFound = DIAHelpers::integrateWindow(spectrum, mz, im, intensity, mz_range, im_range, dia_centroided_); double ppmdiff = Math::getPPMAbs(mz, b_ion_mz); if (signalFound && ppmdiff < dia_byseries_ppm_diff_ && intensity > dia_byseries_intensity_min_) { @@ -256,11 +250,9 @@ namespace OpenMS } for (const auto& y_ion_mz : yseries) { - left = y_ion_mz; - right = y_ion_mz; - DIAHelpers::adjustExtractionWindow(right, left, dia_extract_window_, dia_extraction_ppm_); + RangeMZ mz_range = DIAHelpers::createMZRangePPM(y_ion_mz, dia_extract_window_, dia_extraction_ppm_); - bool signalFound = DIAHelpers::integrateWindow(spectrum, left, right, mz, intensity, dia_centroided_); + bool signalFound = DIAHelpers::integrateWindow(spectrum, mz, im, intensity, mz_range, im_range, dia_centroided_); double ppmdiff = Math::getPPMAbs(mz, y_ion_mz); if (signalFound && ppmdiff < dia_byseries_ppm_diff_ && intensity > dia_byseries_intensity_min_) { @@ -269,11 +261,10 @@ namespace OpenMS } } - void DIAScoring::score_with_isotopes(SpectrumPtrType spectrum, const std::vector& transitions, - double& dotprod, double& manhattan) const + void DIAScoring::score_with_isotopes(SpectrumSequence& spectrum, const std::vector& transitions, const RangeMobility& im_range, double& dotprod, double& manhattan) const { OpenMS::DiaPrescore dp(dia_extract_window_, dia_nr_isotopes_, dia_nr_charges_); - dp.score(std::move(spectrum), transitions, dotprod, manhattan); + dp.score(spectrum, transitions, im_range, dotprod, manhattan); } /////////////////////////////////////////////////////////////////////////// @@ -292,8 +283,9 @@ namespace OpenMS } } - void DIAScoring::diaIsotopeScoresSub_(const std::vector& transitions, const SpectrumPtrType& spectrum, + void DIAScoring::diaIsotopeScoresSub_(const std::vector& transitions, const SpectrumSequence& spectrum, std::map& intensities, //relative intensities + const RangeMobility& im_range, double& isotope_corr, double& isotope_overlap) const { @@ -317,11 +309,9 @@ namespace OpenMS double abs_charge = std::fabs(static_cast(putative_fragment_charge)); for (int iso = 0; iso <= dia_nr_isotopes_; ++iso) { - double left = transitions[k].getProductMZ() + iso * C13C12_MASSDIFF_U / abs_charge; - double right = left; - DIAHelpers::adjustExtractionWindow(right, left, dia_extract_window_, dia_extraction_ppm_); - double mz, intensity; - DIAHelpers::integrateWindow(spectrum, left, right, mz, intensity, dia_centroided_); + RangeMZ mz_range = DIAHelpers::createMZRangePPM(transitions[k].getProductMZ() + iso * C13C12_MASSDIFF_U / abs_charge, dia_extract_window_, dia_extraction_ppm_); + double mz, intensity, im; + DIAHelpers::integrateWindow(spectrum, mz, im, intensity, mz_range, im_range, dia_centroided_); isotopes_int.push_back(intensity); } @@ -329,25 +319,23 @@ namespace OpenMS // isotope correlation (forward) and the isotope overlap (backward) scores double score = scoreIsotopePattern_(isotopes_int, transitions[k].getProductMZ(), putative_fragment_charge); isotope_corr += score * rel_intensity; - largePeaksBeforeFirstIsotope_(spectrum, transitions[k].getProductMZ(), isotopes_int[0], nr_occurences, max_ratio); + largePeaksBeforeFirstIsotope_(spectrum, transitions[k].getProductMZ(), isotopes_int[0], nr_occurences, max_ratio, im_range); isotope_overlap += nr_occurences * rel_intensity; } } - void DIAScoring::largePeaksBeforeFirstIsotope_(const SpectrumPtrType& spectrum, double mono_mz, double mono_int, int& nr_occurences, double& max_ratio) const + void DIAScoring::largePeaksBeforeFirstIsotope_(const SpectrumSequence& spectrum, double mono_mz, double mono_int, int& nr_occurences, double& max_ratio, const RangeMobility& im_range) const { - double mz, intensity; + double mz, intensity, im; nr_occurences = 0; max_ratio = 0.0; for (int ch = 1; ch <= dia_nr_charges_; ++ch) { - double center = mono_mz - C13C12_MASSDIFF_U / (double) ch; - double left = center; - double right = center; - DIAHelpers::adjustExtractionWindow(right, left, dia_extract_window_, dia_extraction_ppm_); - bool signalFound = DIAHelpers::integrateWindow(spectrum, left, right, mz, intensity, dia_centroided_); + double center = mono_mz - C13C12_MASSDIFF_U / (double) ch; + RangeMZ mz_range = DIAHelpers::createMZRangePPM(center, dia_extract_window_, dia_extraction_ppm_); + bool signalFound = DIAHelpers::integrateWindow(spectrum, mz, im, intensity, mz_range, im_range, dia_centroided_); // Continue if no signal was found - we therefore don't make a statement // about the mass difference if no signal is present. if (!signalFound) diff --git a/src/openms/source/ANALYSIS/OPENSWATH/IonMobilityScoring.cpp b/src/openms/source/ANALYSIS/OPENSWATH/IonMobilityScoring.cpp index c3b0fb819eb..53530c82b48 100644 --- a/src/openms/source/ANALYSIS/OPENSWATH/IonMobilityScoring.cpp +++ b/src/openms/source/ANALYSIS/OPENSWATH/IonMobilityScoring.cpp @@ -26,26 +26,13 @@ namespace OpenMS { - - struct MobilityPeak - { - double im; - double intensity; - MobilityPeak (); - MobilityPeak (double im_, double int_) : - im(im_), - intensity(int_) - {} - }; - typedef std::vector< MobilityPeak > IonMobilogram; - - std::vector computeGrid(const std::vector< IonMobilogram >& mobilograms, double eps) + std::vector IonMobilityScoring::computeGrid_(const std::vector< IonMobilogram >& mobilograms, double eps) { // Extract all ion mobility values across all transitions and produce a // grid of all permitted ion mobility values std::vector im_grid; std::vector< double > mobilityValues; - for (const auto & im_profile : mobilograms) + for (const auto & im_profile : mobilograms) { mobilityValues.reserve(mobilityValues.size() + im_profile.size()); for (const auto & k : im_profile) mobilityValues.push_back(k.im); @@ -55,13 +42,13 @@ namespace OpenMS std::sort(mobilityValues.begin(), mobilityValues.end()); // Reduce mobility values to grid (consider equal if closer than eps) - // + // // In some cases there are not enough datapoints available (one of the // transitions has no datapoints) if (!mobilityValues.empty()) { im_grid.push_back( mobilityValues[0] ); - for (Size k = 1; k < mobilityValues.size(); k++) + for (Size k = 1; k < mobilityValues.size(); k++) { double diff = fabs(mobilityValues[k] - mobilityValues[k-1]); if (diff > eps) @@ -89,7 +76,7 @@ namespace OpenMS @param max_peak_idx The grid position of the maximum */ - void alignToGrid(const IonMobilogram& profile, + void IonMobilityScoring::alignToGrid_(const IonMobilogram& profile, const std::vector& im_grid, std::vector< double >& al_int_values, std::vector< double >& al_im_values, @@ -133,28 +120,15 @@ namespace OpenMS } } - /** - @brief Integrate intensity in an ion mobility spectrum from start to end - - This function will integrate the intensity in a spectrum between mz_start - and mz_end, returning the total intensity and an intensity-weighted drift - time value. - - This function also returns the full ion mobility profile in "res". - - @note If there is no signal, mz will be set to -1 and intensity to 0 - */ - void integrateDriftSpectrum(const OpenSwath::SpectrumPtr& spectrum, - double mz_start, - double mz_end, + // compute ion mobilogram as well as im weighted average. This is based off of integrateWindows() in DIAHelper.cpp + void IonMobilityScoring::computeIonMobilogram(const SpectrumSequence& spectra, + const RangeMZ& mz_range, + const RangeMobility& im_range, double & im, double & intensity, - IonMobilogram& res, - double eps, - double drift_start, - double drift_end) + IonMobilogram& res, + double eps) { - OPENMS_PRECONDITION(spectrum->getDriftTimeArray() != nullptr, "Cannot filter by drift time if no drift time is available."); // rounding multiplier for the ion mobility value // TODO: how to improve this -- will work up to 42949.67296 @@ -164,31 +138,43 @@ namespace OpenMS // same spot in the ion mobilogram (they are not sorted by ion mobility in // the input data), therefore create a map to map to bins. std::map< int, double> im_chrom; - auto mz_arr_end = spectrum->getMZArray()->data.end(); - auto int_it = spectrum->getIntensityArray()->data.begin(); - auto im_it = spectrum->getDriftTimeArray()->data.begin(); - - // this assumes that the spectra are sorted! - auto mz_it = std::lower_bound(spectrum->getMZArray()->data.begin(), mz_arr_end, mz_start); - auto mz_it_end = std::lower_bound(mz_it, mz_arr_end, mz_end); - // also advance intensity and ion mobility iterator now - auto iterator_pos = std::distance(spectrum->getMZArray()->data.begin(), mz_it); - std::advance(int_it, iterator_pos); - std::advance(im_it, iterator_pos); - - // Iterate from mz start to end, only storing ion mobility values that are in the range - for (; mz_it != mz_it_end; ++mz_it, ++int_it, ++im_it) + for (const auto& spectrum : spectra) { - if ( *im_it >= drift_start && *im_it <= drift_end) + OPENMS_PRECONDITION(spectrum->getDriftTimeArray() != nullptr, "Cannot filter by drift time if no drift time is available."); + OPENMS_PRECONDITION(spectrum->getMZArray()->data.size() == spectrum->getIntensityArray()->data.size(), "MZ and Intensity array need to have the same length."); + OPENMS_PRECONDITION(spectrum->getMZArray()->data.size() == spectrum->getDriftTimeArray()->data.size(), "MZ and Drift Time array need to have the same length."); + + auto mz_arr_end = spectrum->getMZArray()->data.end(); + auto int_it = spectrum->getIntensityArray()->data.begin(); + auto im_it = spectrum->getDriftTimeArray()->data.begin(); + + // this assumes that the spectra are sorted! + auto mz_it = std::lower_bound(spectrum->getMZArray()->data.begin(), mz_arr_end, mz_range.getMin()); + // auto mz_it_end = std::lower_bound(mz_it, mz_arr_end, mz_end); + + // also advance intensity and ion mobility iterator now + auto iterator_pos = std::distance(spectrum->getMZArray()->data.begin(), mz_it); + std::advance(int_it, iterator_pos); + std::advance(im_it, iterator_pos); + + // Start iteration from mz start, end iteration when mz value is larger than mz_end, only store only storing ion mobility values that are in the range + double mz_end = mz_range.getMax(); + while ( ( *mz_it < mz_end ) && (mz_it < mz_arr_end) ) { - // std::cout << "IM " << *im_it << " mz " << *mz_it << " int " << *int_it << std::endl; - im_chrom[ int((*im_it)*IM_IDX_MULT) ] += *int_it; - intensity += (*int_it); - im += (*int_it) * (*im_it); + if (im_range.contains(*im_it)) + { + intensity += (*int_it); + im += (*int_it) * (*im_it); + im_chrom[ int((*im_it)*IM_IDX_MULT) ] += *int_it; + } + ++mz_it; + ++int_it; + ++im_it; } } + // compute the weighted average ion mobility if (intensity > 0.) { im /= intensity; @@ -212,29 +198,41 @@ namespace OpenMS /// Destructor IonMobilityScoring::~IonMobilityScoring() = default; - void IonMobilityScoring::driftScoringMS1Contrast(const OpenSwath::SpectrumPtr& spectrum, const OpenSwath::SpectrumPtr& ms1spectrum, + void IonMobilityScoring::driftScoringMS1Contrast(const SpectrumSequence& spectra, const SpectrumSequence& ms1spectrum, const std::vector & transitions, OpenSwath_Scores & scores, - const double drift_lower, - const double drift_upper, + RangeMobility im_range, const double dia_extract_window_, const bool dia_extraction_ppm_, const double drift_extra) { - OPENMS_PRECONDITION(spectrum != nullptr, "Spectrum cannot be null"); + OPENMS_PRECONDITION(!spectra.empty(), "Spectra cannot be empty") + OPENMS_PRECONDITION(!ms1spectrum.empty(), "MS1 spectrum cannot be empty") OPENMS_PRECONDITION(!transitions.empty(), "Need at least one transition"); - if (ms1spectrum->getDriftTimeArray() == nullptr) + //TODO not sure what error format is best + for (const auto& s:spectra) { - OPENMS_LOG_DEBUG << " ERROR: Drift time is missing in ion mobility spectrum!" << std::endl; - return; + if (s->getDriftTimeArray() == nullptr) + { + OPENMS_LOG_DEBUG << " ERROR: Drift time is missing in ion mobility spectrum!" << std::endl; + return; + } + } + + for (const auto& s:ms1spectrum) + { + if (s->getDriftTimeArray() == nullptr) + { + OPENMS_LOG_DEBUG << " ERROR: Drift time is missing in MS1 ion mobility spectrum!" << std::endl; + return; + } } double eps = 1e-5; // eps for two grid cells to be considered equal - double drift_width = fabs(drift_upper - drift_lower); - double drift_lower_used = drift_lower - drift_width * drift_extra; - double drift_upper_used = drift_upper + drift_width * drift_extra; + // extend IM range by drift_extra + im_range.scaleBy(drift_extra * 2. + 1); // multiple by 2 because want drift extra to be extended by that amount on either side std::vector< IonMobilogram > mobilograms; @@ -245,43 +243,42 @@ namespace OpenMS IonMobilogram res; const TransitionType transition = transitions[k]; // Calculate the difference of the theoretical ion mobility and the actually measured ion mobility - double left(transition.getProductMZ()), right(transition.getProductMZ()); - DIAHelpers::adjustExtractionWindow(right, left, dia_extract_window_, dia_extraction_ppm_); + RangeMZ mz_range = DIAHelpers::createMZRangePPM(transition.getProductMZ(), dia_extract_window_, dia_extraction_ppm_); - integrateDriftSpectrum(spectrum, left, right, im, intensity, res, eps, drift_lower_used, drift_upper_used); + computeIonMobilogram(spectra, mz_range, im_range, im, intensity, res, eps); mobilograms.push_back( std::move(res) ); } // Step 2: MS1 extraction double im(0), intensity(0); IonMobilogram ms1_profile; - double left(transitions[0].getPrecursorMZ()), right(transitions[0].getPrecursorMZ()); - DIAHelpers::adjustExtractionWindow(right, left, dia_extract_window_, dia_extraction_ppm_); - integrateDriftSpectrum(ms1spectrum, left, right, im, intensity, ms1_profile, eps, drift_lower_used, drift_upper_used); // TODO: aggregate over isotopes + RangeMZ mz_range = DIAHelpers::createMZRangePPM(transitions[0].getPrecursorMZ(), dia_extract_window_, dia_extraction_ppm_); + + computeIonMobilogram(ms1spectrum, mz_range, im_range, im, intensity, ms1_profile, eps); // TODO: aggregate over isotopes mobilograms.push_back(ms1_profile); - std::vector im_grid = computeGrid(mobilograms, eps); // ensure grid is based on all profiles! + std::vector im_grid = computeGrid_(mobilograms, eps); // ensure grid is based on all profiles! mobilograms.pop_back(); // Step 3: Align the IonMobilogram vectors to the grid std::vector< std::vector< double > > aligned_mobilograms; - for (const auto & mobilogram : mobilograms) + for (const auto & mobilogram : mobilograms) { std::vector< double > arrInt, arrIM; Size max_peak_idx = 0; - alignToGrid(mobilogram, im_grid, arrInt, arrIM, eps, max_peak_idx); + alignToGrid_(mobilogram, im_grid, arrInt, arrIM, eps, max_peak_idx); aligned_mobilograms.push_back(arrInt); } std::vector< double > ms1_int_values, ms1_im_values; Size max_peak_idx = 0; - alignToGrid(ms1_profile, im_grid, ms1_int_values, ms1_im_values, eps, max_peak_idx); + alignToGrid_(ms1_profile, im_grid, ms1_int_values, ms1_im_values, eps, max_peak_idx); // Step 4: MS1 contrast scores { OpenSwath::MRMScoring mrmscore_; mrmscore_.initializeXCorrPrecursorContrastMatrix({ms1_int_values}, aligned_mobilograms); - OPENMS_LOG_DEBUG << "all-all: Contrast Scores : coelution precursor : " << mrmscore_.calcXcorrPrecursorContrastCoelutionScore() << " / shape precursor " << + OPENMS_LOG_DEBUG << "all-all: Contrast Scores : coelution precursor : " << mrmscore_.calcXcorrPrecursorContrastCoelutionScore() << " / shape precursor " << mrmscore_.calcXcorrPrecursorContrastShapeScore() << std::endl; scores.im_ms1_contrast_coelution = mrmscore_.calcXcorrPrecursorContrastCoelutionScore(); scores.im_ms1_contrast_shape = mrmscore_.calcXcorrPrecursorContrastShapeScore(); @@ -301,7 +298,7 @@ namespace OpenMS OpenSwath::MRMScoring mrmscore_; // horribly broken: provides vector of length 1, but expects at least length 2 in calcXcorrPrecursorContrastCoelutionScore() mrmscore_.initializeXCorrPrecursorContrastMatrix({ms1_int_values}, {fragment_values}); - OPENMS_LOG_DEBUG << "Contrast Scores : coelution precursor : " << mrmscore_.calcXcorrPrecursorContrastSumFragCoelutionScore() << " / shape precursor " << + OPENMS_LOG_DEBUG << "Contrast Scores : coelution precursor : " << mrmscore_.calcXcorrPrecursorContrastSumFragCoelutionScore() << " / shape precursor " << mrmscore_.calcXcorrPrecursorContrastSumFragShapeScore() << std::endl; // in order to prevent assertion error call calcXcorrPrecursorContrastSumFragCoelutionScore, same as calcXcorrPrecursorContrastCoelutionScore() however different assertion @@ -309,37 +306,35 @@ namespace OpenMS // in order to prevent assertion error call calcXcorrPrecursorContrastSumFragShapeScore(), same as calcXcorrPrecursorContrastShapeScore() however different assertion. scores.im_ms1_sum_contrast_shape = mrmscore_.calcXcorrPrecursorContrastSumFragShapeScore(); - } - void IonMobilityScoring::driftScoringMS1(const OpenSwath::SpectrumPtr& spectrum, + void IonMobilityScoring::driftScoringMS1(const SpectrumSequence & spectra, const std::vector & transitions, OpenSwath_Scores & scores, - const double drift_lower, - const double drift_upper, const double drift_target, + RangeMobility im_range, const double dia_extract_window_, const bool dia_extraction_ppm_, const bool /* use_spline */, const double drift_extra) { - OPENMS_PRECONDITION(spectrum != nullptr, "Spectrum cannot be null"); + OPENMS_PRECONDITION(!spectra.empty(), "Spectra cannot be empty") OPENMS_PRECONDITION(!transitions.empty(), "Need at least one transition"); - if (spectrum->getDriftTimeArray() == nullptr) - { - OPENMS_LOG_DEBUG << " ERROR: Drift time is missing in ion mobility spectrum!" << std::endl; - return; + for (auto s:spectra){ + if (s->getDriftTimeArray() == nullptr) + { + OPENMS_LOG_DEBUG << " ERROR: Drift time is missing in ion mobility spectrum!" << std::endl; + return; + } } - double drift_width = fabs(drift_upper - drift_lower); - double drift_lower_used = drift_lower - drift_width * drift_extra; - double drift_upper_used = drift_upper + drift_width * drift_extra; + im_range.scaleBy(drift_extra * 2. + 1); // multiple by 2 because want drift extra to be extended by that amount on either side - double im(0), intensity(0); - double left(transitions[0].getPrecursorMZ()), right(transitions[0].getPrecursorMZ()); - DIAHelpers::adjustExtractionWindow(right, left, dia_extract_window_, dia_extraction_ppm_); - DIAHelpers::integrateDriftSpectrum(spectrum, left, right, im, intensity, drift_lower_used, drift_upper_used); + double im(0), intensity(0), mz(0); + RangeMZ mz_range = DIAHelpers::createMZRangePPM(transitions[0].getPrecursorMZ(), dia_extract_window_, dia_extraction_ppm_); + + DIAHelpers::integrateWindow(spectra, mz, im, intensity, mz_range, im_range); // Record the measured ion mobility scores.im_ms1_drift = im; @@ -349,30 +344,29 @@ namespace OpenMS scores.im_ms1_delta = drift_target - im; } - void IonMobilityScoring::driftScoring(const OpenSwath::SpectrumPtr& spectrum, + void IonMobilityScoring::driftScoring(const SpectrumSequence& spectra, const std::vector & transitions, OpenSwath_Scores & scores, - const double drift_lower, - const double drift_upper, const double drift_target, + RangeMobility im_range, const double dia_extract_window_, const bool dia_extraction_ppm_, const bool /* use_spline */, const double drift_extra) { - OPENMS_PRECONDITION(spectrum != nullptr, "Spectrum cannot be null"); - - if (spectrum->getDriftTimeArray() == nullptr) + OPENMS_PRECONDITION(!spectra.empty(), "Spectra cannot be empty"); + for (auto s:spectra) { - OPENMS_LOG_DEBUG << " ERROR: Drift time is missing in ion mobility spectrum!" << std::endl; - return; + if (s->getDriftTimeArray() == nullptr) + { + OPENMS_LOG_DEBUG << " ERROR: Drift time is missing in ion mobility spectrum!" << std::endl; + return; + } } double eps = 1e-5; // eps for two grid cells to be considered equal - double drift_width = fabs(drift_upper - drift_lower); - double drift_lower_used = drift_lower - drift_width * drift_extra; - double drift_upper_used = drift_upper + drift_width * drift_extra; + im_range.scaleBy(drift_extra * 2. + 1); // multiple by 2 because want drift extra to be extended by that amount on either side double delta_drift = 0; double delta_drift_abs = 0; @@ -391,9 +385,11 @@ namespace OpenMS double im(0), intensity(0); // Calculate the difference of the theoretical ion mobility and the actually measured ion mobility - double left(transition.getProductMZ()), right(transition.getProductMZ()); - DIAHelpers::adjustExtractionWindow(right, left, dia_extract_window_, dia_extraction_ppm_); - integrateDriftSpectrum(spectrum, left, right, im, intensity, res, eps, drift_lower_used, drift_upper_used); + RangeMZ mz_range = DIAHelpers::createMZRangePPM(transition.getProductMZ(), dia_extract_window_, dia_extraction_ppm_); + + //double left(transition.getProductMZ()), right(transition.getProductMZ()); + //DIAHelpers::adjustExtractionWindow(right, left, dia_extract_window_, dia_extraction_ppm_); + computeIonMobilogram(spectra, mz_range, im_range, im, intensity, res, eps); mobilograms.push_back(res); // TODO what do to about those that have no signal ? @@ -434,13 +430,13 @@ namespace OpenMS scores.im_drift_weighted = computed_im_weighted; // Step 2: Align the IonMobilogram vectors to the grid - std::vector im_grid = computeGrid(mobilograms, eps); + std::vector im_grid = computeGrid_(mobilograms, eps); std::vector< std::vector< double > > aligned_mobilograms; - for (const auto & mobilogram : mobilograms) + for (const auto & mobilogram : mobilograms) { std::vector< double > arr_int, arr_IM; Size max_peak_idx = 0; - alignToGrid(mobilogram, im_grid, arr_int, arr_IM, eps, max_peak_idx); + alignToGrid_(mobilogram, im_grid, arr_int, arr_IM, eps, max_peak_idx); if (!arr_int.empty()) aligned_mobilograms.push_back(arr_int); } @@ -451,6 +447,8 @@ namespace OpenMS scores.im_xcorr_shape_score = std::numeric_limits::quiet_NaN(); return; } + + OpenSwath::MRMScoring mrmscore_; mrmscore_.initializeXCorrMatrix(aligned_mobilograms); @@ -460,6 +458,4 @@ namespace OpenMS scores.im_xcorr_coelution_score = xcorr_coelution_score; scores.im_xcorr_shape_score = xcorr_shape_score; } - } - diff --git a/src/openms/source/ANALYSIS/OPENSWATH/MRMFeatureFinderScoring.cpp b/src/openms/source/ANALYSIS/OPENSWATH/MRMFeatureFinderScoring.cpp index 78fcf0394f8..fa936706ff5 100644 --- a/src/openms/source/ANALYSIS/OPENSWATH/MRMFeatureFinderScoring.cpp +++ b/src/openms/source/ANALYSIS/OPENSWATH/MRMFeatureFinderScoring.cpp @@ -84,6 +84,7 @@ namespace OpenMS defaults_.setMinFloat("im_extra_drift", 0.0); defaults_.setValue("strict", "true", "Whether to error (true) or skip (false) if a transition in a transition group does not have a corresponding chromatogram.", {"advanced"}); defaults_.setValidStrings("strict", {"true","false"}); + defaults_.setValue("use_ms1_ion_mobility", "true", "Performs ion mobility extraction in MS1. Set to false if MS1 spectra do not contain ion mobility", {"advanced"}); defaults_.insert("TransitionGroupPicker:", MRMTransitionGroupPicker().getDefaults()); @@ -168,7 +169,7 @@ namespace OpenMS void MRMFeatureFinderScoring::pickExperiment(const OpenSwath::SpectrumAccessPtr& input, FeatureMap& output, const OpenSwath::LightTargetedExperiment& transition_exp, - const TransformationDescription& trafo, + const TransformationDescription& trafo, const std::vector& swath_maps, TransitionGroupMapType& transition_group_map) { @@ -305,22 +306,17 @@ namespace OpenMS { MRMFeature idmrmfeature = trgr_ident.getFeaturesMuteable()[feature_idx]; OpenSwath::IMRMFeature* idimrmfeature; - idimrmfeature = new MRMFeatureOpenMS(idmrmfeature); + idimrmfeature = new MRMFeatureOpenMS(idmrmfeature); // get drift time upper/lower offset (this assumes that all chromatograms // are derived from the same precursor with the same drift time) - double drift_lower(0), drift_upper(0); - if (!trgr_ident.getChromatograms().empty()) + RangeMobility im_range; + + if ( (!trgr_ident.getChromatograms().empty()) || (!trgr_ident.getPrecursorChromatograms().empty()) ) { auto & prec = trgr_ident.getChromatograms()[0].getPrecursor(); - drift_lower = prec.getDriftTime() - prec.getDriftTimeWindowLowerOffset(); - drift_upper = prec.getDriftTime() + prec.getDriftTimeWindowUpperOffset(); - } - else if (!trgr_ident.getPrecursorChromatograms().empty()) - { - auto & prec = trgr_ident.getPrecursorChromatograms()[0].getPrecursor(); - drift_lower = prec.getDriftTime() - prec.getDriftTimeWindowLowerOffset(); - drift_upper = prec.getDriftTime() + prec.getDriftTimeWindowUpperOffset(); + im_range.setMin(prec.getDriftTime()); // sets the minimum and maximum + im_range.minSpanIfSingular(prec.getDriftTimeWindowLowerOffset()); } std::vector native_ids_identification; @@ -331,7 +327,7 @@ namespace OpenMS OpenSwath::ISignalToNoisePtr snptr(new OpenMS::SignalToNoiseOpenMS< MSChromatogram >( trgr_ident.getChromatogram(trgr_ident.getTransitions()[i].getNativeID()), sn_win_len_, sn_bin_count_, write_log_messages_)); - if ( (snptr->getValueAtRT(idmrmfeature.getRT()) > uis_threshold_sn_) + if ( (snptr->getValueAtRT(idmrmfeature.getRT()) > uis_threshold_sn_) && (idmrmfeature.getFeature(trgr_ident.getTransitions()[i].getNativeID()).getIntensity() > uis_threshold_peak_area_)) { signal_noise_estimators_identification.push_back(snptr); @@ -343,7 +339,7 @@ namespace OpenMS if (!native_ids_identification.empty()) { scorer.calculateChromatographicIdScores(idimrmfeature, - native_ids_identification, + native_ids_identification, native_ids_detection, signal_noise_estimators_identification, idscores); @@ -430,9 +426,9 @@ namespace OpenMS { OpenSwath_Scores tmp_scores; - scorer.calculateDIAIdScores(idimrmfeature, + scorer.calculateDIAIdScores(idimrmfeature, trgr_ident.getTransition(native_ids_identification[i]), - swath_maps, diascoring_, tmp_scores, drift_lower, drift_upper); + swath_maps, im_range, diascoring_, tmp_scores); ind_isotope_correlation.push_back(tmp_scores.isotope_correlation); ind_isotope_overlap.push_back(tmp_scores.isotope_overlap); @@ -448,9 +444,9 @@ namespace OpenMS } void MRMFeatureFinderScoring::scorePeakgroups(MRMTransitionGroupType& transition_group, - const TransformationDescription& trafo, + const TransformationDescription& trafo, const std::vector& swath_maps, - FeatureMap& output, + FeatureMap& output, bool ms1only) const { if (PeptideRefMap_.empty()) @@ -476,20 +472,29 @@ namespace OpenMS // get drift time upper/lower offset (this assumes that all chromatograms // are derived from the same precursor with the same drift time) - double drift_lower(0), drift_upper(0), drift_target(0); - if (!transition_group_detection.getChromatograms().empty()) + RangeMobility im_range; + double drift_target(0); + + auto setDriftTarget = [](auto& prec){ + double lower_bound = prec.getDriftTime() - prec.getDriftTimeWindowLowerOffset(); + double upper_bound = prec.getDriftTime() + prec.getDriftTimeWindowUpperOffset(); + return RangeMobility(lower_bound, upper_bound); + }; + + if ( !transition_group_detection.getChromatograms().empty() ) { auto & prec = transition_group_detection.getChromatograms()[0].getPrecursor(); - drift_lower = prec.getDriftTime() - prec.getDriftTimeWindowLowerOffset(); - drift_upper = prec.getDriftTime() + prec.getDriftTimeWindowUpperOffset(); drift_target = prec.getDriftTime(); + + if (drift_target > 0) im_range = setDriftTarget(prec); } - else if (!transition_group_detection.getPrecursorChromatograms().empty()) + + else if ( !transition_group_detection.getPrecursorChromatograms().empty() ) { auto & prec = transition_group_detection.getPrecursorChromatograms()[0].getPrecursor(); - drift_lower = prec.getDriftTime() - prec.getDriftTimeWindowLowerOffset(); - drift_upper = prec.getDriftTime() + prec.getDriftTimeWindowUpperOffset(); drift_target = prec.getDriftTime(); + + if (drift_target > 0) im_range = setDriftTarget(prec); } // currently we cannot do much about the log messages and they mostly occur in decoy transition signals @@ -524,7 +529,8 @@ namespace OpenMS spacing_for_spectra_resampling_, im_extra_drift_, su_, - spectrum_addition_method_); + spectrum_addition_method_, + use_ms1_ion_mobility_); ProteaseDigestion pd; pd.setEnzyme("Trypsin"); @@ -544,7 +550,7 @@ namespace OpenMS OPENMS_LOG_DEBUG << "Scoring feature " << (mrmfeature) << " == " << mrmfeature.getMetaValue("PeptideRef") << " [ expected RT " << PeptideRefMap_.at(mrmfeature.getMetaValue("PeptideRef"))->rt << " / " << expected_rt << " ]" << - " with " << transition_group_detection.size() << " transitions and " << + " with " << transition_group_detection.size() << " transitions and " << transition_group_detection.getChromatograms().size() << " chromatograms" << std::endl; int group_size = boost::numeric_cast(transition_group_detection.size()); @@ -552,7 +558,7 @@ namespace OpenMS { delete imrmfeature; // free resources before continuing throw Exception::IllegalArgument(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, - "Error: Transition group " + transition_group_detection.getTransitionGroupID() + + "Error: Transition group " + transition_group_detection.getTransitionGroupID() + " has no chromatograms."); } bool swath_present = (!swath_maps.empty() && swath_maps[0].sptr->getNrSpectra() > 0); @@ -574,15 +580,15 @@ namespace OpenMS scores.sn_ratio = mrmscore_.calcSNScore(imrmfeature, ms1_signal_noise_estimators); // everything below S/N 1 can be set to zero (and the log safely applied) if (scores.sn_ratio < 1) - { + { scores.log_sn_score = 0; } else - { + { scores.log_sn_score = std::log(scores.sn_ratio); } - if (su_.use_sn_score_) - { + if (su_.use_sn_score_) + { mrmfeature.addScore("sn_ratio", scores.sn_ratio); mrmfeature.addScore("var_log_sn_score", scores.log_sn_score); // compute subfeature log-SN values @@ -614,10 +620,10 @@ namespace OpenMS mrmfeature.addScore("var_norm_rt_score", scores.norm_rt_score); } - // full spectra scores + // full spectra scores if (ms1_map_ && ms1_map_->getNrSpectra() > 0 && mrmfeature.getMZ() > 0) { - scorer.calculatePrecursorDIAScores(ms1_map_, diascoring_, precursor_mz, imrmfeature->getRT(), *pep, scores, drift_lower, drift_upper); + scorer.calculatePrecursorDIAScores(ms1_map_, diascoring_, precursor_mz, imrmfeature->getRT(), *pep, im_range, scores); } if (su_.use_ms1_fullscan) { @@ -676,7 +682,7 @@ namespace OpenMS scorer.calculateDIAScores(imrmfeature, transition_group_detection.getTransitions(), swath_maps, ms1_map_, diascoring_, *pep, scores, masserror_ppm, - drift_lower, drift_upper, drift_target); + drift_target, im_range); mrmfeature.setMetaValue("masserror_ppm", masserror_ppm); } if (sonar_present && su_.use_sonar_scores) @@ -930,7 +936,7 @@ namespace OpenMS mrmfeature.setMetaValue("missedCleavages", pd.peptideCount(pep_hit_.getSequence()) - 1); } - // set protein accession numbers + // set protein accession numbers for (Size k = 0; k < pep->protein_refs.size(); k++) { PeptideEvidence pe; @@ -1030,8 +1036,10 @@ namespace OpenMS sonarscoring_.setParameters(p); diascoring_.setParameters(param_.copy("DIAScoring:", true)); + emgscoring_.setFitterParam(param_.copy("EMGScoring:", true)); strict_ = (bool)param_.getValue("strict").toBool(); + use_ms1_ion_mobility_ = (bool)param_.getValue("use_ms1_ion_mobility").toBool(); su_.use_coelution_score_ = param_.getValue("Scores:use_coelution_score").toBool(); su_.use_shape_score_ = param_.getValue("Scores:use_shape_score").toBool(); diff --git a/src/openms/source/ANALYSIS/OPENSWATH/OpenSwathScoring.cpp b/src/openms/source/ANALYSIS/OPENSWATH/OpenSwathScoring.cpp index 66b8a6eb48b..82f4f7634f2 100644 --- a/src/openms/source/ANALYSIS/OPENSWATH/OpenSwathScoring.cpp +++ b/src/openms/source/ANALYSIS/OPENSWATH/OpenSwathScoring.cpp @@ -24,50 +24,6 @@ #include #include - -namespace OpenMS -{ - - void sortSpectrumByMZ(OpenSwath::Spectrum& spec) - { - //sort index list - std::vector > sorted_indices; - sorted_indices.reserve(spec.getMZArray()->data.size()); - auto mz_it = spec.getMZArray()->data.begin(); - for (Size i = 0; i < spec.getMZArray()->data.size(); ++i) - { - sorted_indices.emplace_back(*mz_it, i); - ++mz_it; - } - std::stable_sort(sorted_indices.begin(), sorted_indices.end()); - - // extract list of indices - std::vector select_indices; - select_indices.reserve(sorted_indices.size()); - for (const auto& sidx : sorted_indices) - { - select_indices.push_back(sidx.second); - } - - for (auto& da : spec.getDataArrays() ) - { - if (da->data.empty()) continue; - OpenSwath::BinaryDataArrayPtr tmp(new OpenSwath::BinaryDataArray); - tmp->description = da->description; - tmp->data.reserve(select_indices.size()); - for (Size i = 0; i < select_indices.size(); ++i) - { - tmp->data.push_back( da->data[ select_indices[i] ] ); - } - da = tmp; - } - - OPENMS_POSTCONDITION( std::adjacent_find(spec.getMZArray()->data.begin(), - spec.getMZArray()->data.end(), std::greater()) == spec.getMZArray()->data.end(), - "Postcondition violated: m/z vector needs to be sorted!" ) - } -} - namespace OpenMS { @@ -76,7 +32,7 @@ namespace OpenMS rt_normalization_factor_(1.0), spacing_for_spectra_resampling_(0.005), add_up_spectra_(1), - spectra_addition_method_("simple"), + spectra_addition_method_(SpectrumAdditionMethod::ADDITION), im_drift_extra_pcnt_(0.0) { } @@ -89,14 +45,28 @@ namespace OpenMS double spacing_for_spectra_resampling, const double drift_extra, const OpenSwath_Scores_Usage & su, - const std::string& spectrum_addition_method) + const std::string& spectrum_addition_method, + bool use_ms1_ion_mobility) { this->rt_normalization_factor_ = rt_normalization_factor; this->add_up_spectra_ = add_up_spectra; - this->spectra_addition_method_ = spectrum_addition_method; + if (spectrum_addition_method == "simple") + { + this->spectra_addition_method_ = SpectrumAdditionMethod::ADDITION; + } + else if (spectrum_addition_method == "resample") + { + this->spectra_addition_method_ = SpectrumAdditionMethod::RESAMPLE; + } + else + { + throw Exception::InvalidValue(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "spectrum_addition_method must be simple or resample", spectrum_addition_method); + } + this->im_drift_extra_pcnt_ = drift_extra; this->spacing_for_spectra_resampling_ = spacing_for_spectra_resampling; this->su_ = su; + this->use_ms1_ion_mobility_ = use_ms1_ion_mobility; } void OpenSwathScoring::calculateDIAScores(OpenSwath::IMRMFeature* imrmfeature, @@ -107,9 +77,8 @@ namespace OpenMS const CompoundType& compound, OpenSwath_Scores& scores, std::vector& masserror_ppm, - const double drift_lower, - const double drift_upper, - const double drift_target) + const double drift_target,// TODO is this needed + const RangeMobility& im_range) { OPENMS_PRECONDITION(imrmfeature != nullptr, "Feature to be scored cannot be null"); OPENMS_PRECONDITION(transitions.size() > 0, "There needs to be at least one transition."); @@ -138,33 +107,31 @@ namespace OpenMS getNormalized_library_intensities_(transitions, normalized_library_intensity); // find spectrum that is closest to the apex of the peak using binary search - OpenSwath::SpectrumPtr spectrum = fetchSpectrumSwath(used_swath_maps, imrmfeature->getRT(), add_up_spectra_, drift_lower, drift_upper); + std::vector spectra = fetchSpectrumSwath(used_swath_maps, imrmfeature->getRT(), add_up_spectra_, im_range); - // calculate drift extraction width for current spectrum (with some extra for cross-correlation) - double drift_width = fabs(drift_upper - drift_lower); - double drift_lower_used = drift_lower - drift_width * im_drift_extra_pcnt_; - double drift_upper_used = drift_upper + drift_width * im_drift_extra_pcnt_; + // set the DIA parameters + // TODO Cache these parameters + double dia_extract_window_ = (double)diascoring.getParameters().getValue("dia_extraction_window"); + bool dia_extraction_ppm_ = diascoring.getParameters().getValue("dia_extraction_unit") == "ppm"; // score drift time dimension - if (drift_upper > 0 && su_.use_im_scores) + if ( su_.use_im_scores) { - double dia_extract_window_ = (double)diascoring.getParameters().getValue("dia_extraction_window"); - bool dia_extraction_ppm_ = diascoring.getParameters().getValue("dia_extraction_unit") == "ppm"; - auto drift_spectrum = fetchSpectrumSwath(used_swath_maps, imrmfeature->getRT(), add_up_spectra_, drift_lower_used, drift_upper_used); - IonMobilityScoring::driftScoring(drift_spectrum, transitions, scores, - drift_lower, drift_upper, drift_target, + IonMobilityScoring::driftScoring(spectra, transitions, scores, + drift_target, im_range, dia_extract_window_, dia_extraction_ppm_, false, im_drift_extra_pcnt_); } + // Mass deviation score - diascoring.dia_massdiff_score(transitions, spectrum, normalized_library_intensity, scores.massdev_score, scores.weighted_massdev_score, masserror_ppm); + diascoring.dia_massdiff_score(transitions, spectra, normalized_library_intensity, im_range, scores.massdev_score, scores.weighted_massdev_score, masserror_ppm); //TODO this score and the next, both rely on the CoarseIsotope of the PeptideAveragine. Maybe we could // DIA dotproduct and manhattan score based on library intensity and sum formula if present if (su_.use_ms2_isotope_scores) { - diascoring.score_with_isotopes(spectrum, transitions, scores.dotprod_score_dia, scores.manhatt_score_dia); + diascoring.score_with_isotopes(spectra, transitions, im_range, scores.dotprod_score_dia, scores.manhatt_score_dia); // Isotope correlation / overlap score: Is this peak part of an // isotopic pattern or is it the monoisotopic peak in an isotopic @@ -173,7 +140,7 @@ namespace OpenMS // not optimal for metabolites - but better than nothing, given that for // most fragments we don't really know their composition diascoring - .dia_isotope_scores(transitions, spectrum, imrmfeature, scores.isotope_correlation, scores.isotope_overlap); + .dia_isotope_scores(transitions, spectra, imrmfeature, im_range, scores.isotope_correlation, scores.isotope_overlap); } // Peptide-specific scores (only useful, when product transitions are REAL fragments, e.g. not in FFID) @@ -184,51 +151,69 @@ namespace OpenMS OpenMS::AASequence aas; int by_charge_state = 1; // for which charge states should we check b/y series OpenSwathDataAccessHelper::convertPeptideToAASequence(compound, aas); - diascoring.dia_by_ion_score(spectrum, aas, by_charge_state, scores.bseries_score, scores.yseries_score); + diascoring.dia_by_ion_score(spectra, aas, by_charge_state, im_range, scores.bseries_score, scores.yseries_score); + } + + + RangeMobility im_range_ms1; + if (use_ms1_ion_mobility_) + { + im_range_ms1 = im_range; + } + else // do not extract across IM in MS1 + { + im_range_ms1 = RangeMobility(); } - if (ms1_map && ms1_map->getNrSpectra() > 0) + if (ms1_map && ms1_map->getNrSpectra() > 0) { double precursor_mz = transitions[0].precursor_mz; double rt = imrmfeature->getRT(); - calculatePrecursorDIAScores(ms1_map, diascoring, precursor_mz, rt, compound, scores, drift_lower, drift_upper); + calculatePrecursorDIAScores(ms1_map, diascoring, precursor_mz, rt, compound, im_range_ms1, scores); + } + - if (drift_upper > 0 && su_.use_im_scores) - { + if ( (ms1_map && ms1_map->getNrSpectra() > 0) && ( su_.use_im_scores) ) // IM MS1 scores + { double dia_extract_window_ = (double)diascoring.getParameters().getValue("dia_extraction_window"); bool dia_extraction_ppm_ = diascoring.getParameters().getValue("dia_extraction_unit") == "ppm"; - IonMobilityScoring::driftScoringMS1( fetchSpectrumSwath(ms1_map, imrmfeature->getRT(), add_up_spectra_, drift_lower_used, drift_upper_used), - transitions, scores, drift_lower, drift_upper, drift_target, dia_extract_window_, dia_extraction_ppm_, false, im_drift_extra_pcnt_); + double rt = imrmfeature->getRT(); - IonMobilityScoring::driftScoringMS1Contrast( - fetchSpectrumSwath(used_swath_maps, imrmfeature->getRT(), add_up_spectra_, drift_lower_used, drift_upper_used), - fetchSpectrumSwath(ms1_map, imrmfeature->getRT(), add_up_spectra_, drift_lower, drift_upper), - transitions, scores, drift_lower, drift_upper, dia_extract_window_, dia_extraction_ppm_, im_drift_extra_pcnt_); - } - } + std::vector ms1_spectrum = fetchSpectrumSwath(ms1_map, rt, add_up_spectra_, im_range_ms1); + IonMobilityScoring::driftScoringMS1(ms1_spectrum, + transitions, scores, drift_target, im_range_ms1, dia_extract_window_, dia_extraction_ppm_, false, im_drift_extra_pcnt_); + IonMobilityScoring::driftScoringMS1Contrast(spectra, ms1_spectrum, + transitions, scores, im_range_ms1, dia_extract_window_, dia_extraction_ppm_, im_drift_extra_pcnt_); + } } - void OpenSwathScoring::calculatePrecursorDIAScores(const OpenSwath::SpectrumAccessPtr& ms1_map, + void OpenSwathScoring::calculatePrecursorDIAScores(const OpenSwath::SpectrumAccessPtr& ms1_map, const OpenMS::DIAScoring & diascoring, - double precursor_mz, - double rt, - const CompoundType& compound, - OpenSwath_Scores & scores, - double drift_lower, double drift_upper) + double precursor_mz, + double rt, + const CompoundType& compound, + RangeMobility im_range, + OpenSwath_Scores & scores) { + // change im_range based on ms1 settings + if (!use_ms1_ion_mobility_) + { + im_range.clear(); + } + // Compute precursor-level scores: // - compute mass difference in ppm // - compute isotopic pattern score if (ms1_map && ms1_map->getNrSpectra() > 0) { - OpenSwath::SpectrumPtr ms1_spectrum = fetchSpectrumSwath(ms1_map, rt, add_up_spectra_, drift_lower, drift_upper); - diascoring.dia_ms1_massdiff_score(precursor_mz, ms1_spectrum, scores.ms1_ppm_score); + std::vector ms1_spectrum = fetchSpectrumSwath(ms1_map, rt, add_up_spectra_, im_range); + diascoring.dia_ms1_massdiff_score(precursor_mz, ms1_spectrum, im_range, scores.ms1_ppm_score); // derive precursor charge state (get from data if possible) int precursor_charge = 1; - if (compound.getChargeState() != 0) + if (compound.getChargeState() != 0) { precursor_charge = compound.getChargeState(); } @@ -237,15 +222,15 @@ namespace OpenMS { if (!compound.sequence.empty()) { - diascoring.dia_ms1_isotope_scores(precursor_mz, ms1_spectrum, scores.ms1_isotope_correlation, + diascoring.dia_ms1_isotope_scores(precursor_mz, ms1_spectrum, im_range, scores.ms1_isotope_correlation, scores.ms1_isotope_overlap, AASequence::fromString(compound.sequence).getFormula(Residue::Full, precursor_charge)); } else { - diascoring.dia_ms1_isotope_scores_averagine(precursor_mz, ms1_spectrum, + diascoring.dia_ms1_isotope_scores_averagine(precursor_mz, ms1_spectrum, precursor_charge, im_range, scores.ms1_isotope_correlation, - scores.ms1_isotope_overlap, precursor_charge); + scores.ms1_isotope_overlap); } } else @@ -258,15 +243,15 @@ namespace OpenMS // it is shadowed by the exact precursor_mz. //TODO check if charges are the same (in case the charge was actually present in the sum_formula?) empf.setCharge(precursor_charge); - diascoring.dia_ms1_isotope_scores(precursor_mz, ms1_spectrum, scores.ms1_isotope_correlation, + diascoring.dia_ms1_isotope_scores(precursor_mz, ms1_spectrum, im_range, scores.ms1_isotope_correlation, scores.ms1_isotope_overlap, empf); } else { - diascoring.dia_ms1_isotope_scores_averagine(precursor_mz, ms1_spectrum, + diascoring.dia_ms1_isotope_scores_averagine(precursor_mz, ms1_spectrum, precursor_charge, im_range, scores.ms1_isotope_correlation, - scores.ms1_isotope_overlap, precursor_charge); + scores.ms1_isotope_overlap); } } } @@ -275,13 +260,18 @@ namespace OpenMS void OpenSwathScoring::calculateDIAIdScores(OpenSwath::IMRMFeature* imrmfeature, const TransitionType & transition, const std::vector& swath_maps, + RangeMobility& im_range, const OpenMS::DIAScoring & diascoring, - OpenSwath_Scores & scores, - double drift_lower, double drift_upper) + OpenSwath_Scores & scores) { OPENMS_PRECONDITION(imrmfeature != nullptr, "Feature to be scored cannot be null"); OPENMS_PRECONDITION(swath_maps.size() > 0, "There needs to be at least one swath map."); + if (!use_ms1_ion_mobility_) + { + im_range.clear(); + } + // Identify corresponding SONAR maps (if more than one map is used) std::vector used_swath_maps; if (swath_maps.size() > 1) @@ -302,7 +292,7 @@ namespace OpenMS } // find spectrum that is closest to the apex of the peak using binary search - OpenSwath::SpectrumPtr spectrum = fetchSpectrumSwath(used_swath_maps, imrmfeature->getRT(), add_up_spectra_, drift_lower, drift_upper); + std::vector spectrum = fetchSpectrumSwath(used_swath_maps, imrmfeature->getRT(), add_up_spectra_, im_range); // If no charge is given, we assume it to be 1 int putative_product_charge = 1; @@ -316,11 +306,12 @@ namespace OpenMS // pattern? diascoring.dia_ms1_isotope_scores_averagine(transition.getProductMZ(), spectrum, + putative_product_charge, + im_range, scores.isotope_correlation, - scores.isotope_overlap, - putative_product_charge); + scores.isotope_overlap); // Mass deviation score - diascoring.dia_ms1_massdiff_score(transition.getProductMZ(), spectrum, scores.massdev_score); + diascoring.dia_ms1_massdiff_score(transition.getProductMZ(), spectrum, im_range, scores.massdev_score); } void OpenSwathScoring::calculateChromatographicScores( @@ -383,7 +374,7 @@ namespace OpenMS scores.sn_ratio = mrmscore_.calcSNScore(imrmfeature, signal_noise_estimators); // everything below S/N 1 can be set to zero (and the log safely applied) if (scores.sn_ratio < 1) - { + { scores.log_sn_score = 0; } else @@ -506,172 +497,81 @@ namespace OpenMS OpenSwath::Scoring::normalize_sum(&normalized_library_intensity[0], boost::numeric_cast(normalized_library_intensity.size())); } - OpenSwath::SpectrumPtr OpenSwathScoring::fetchSpectrumSwath(OpenSwath::SpectrumAccessPtr swath_map, - double RT, int nr_spectra_to_add, const double drift_lower, const double drift_upper) + SpectrumSequence OpenSwathScoring::fetchSpectrumSwath(OpenSwath::SpectrumAccessPtr swathmap, double RT, int nr_spectra_to_add, const RangeMobility& im_range) { - return getAddedSpectra_(std::move(swath_map), RT, nr_spectra_to_add, drift_lower, drift_upper); - } - OpenSwath::SpectrumPtr OpenSwathScoring::fetchSpectrumSwath(std::vector swath_maps, - double RT, int nr_spectra_to_add, const double drift_lower, const double drift_upper) - { - if (swath_maps.size() == 1) - { - return getAddedSpectra_(swath_maps[0].sptr, RT, nr_spectra_to_add, drift_lower, drift_upper); - } - else + SpectrumSequence all_spectra = swathmap->getMultipleSpectra(RT, nr_spectra_to_add); + if (spectra_addition_method_ == SpectrumAdditionMethod::ADDITION) { - // multiple SWATH maps for a single precursor -> this is SONAR data - std::vector all_spectra; - for (size_t i = 0; i < swath_maps.size(); ++i) - { - OpenSwath::SpectrumPtr spec = getAddedSpectra_(swath_maps[i].sptr, RT, nr_spectra_to_add, drift_lower, drift_upper); - all_spectra.push_back(spec); - } - OpenSwath::SpectrumPtr spectrum_ = SpectrumAddition::addUpSpectra(all_spectra, spacing_for_spectra_resampling_, true); - return spectrum_; + return all_spectra; // return vector, addition is done later } - } - - OpenSwath::SpectrumPtr filterByDrift(const OpenSwath::SpectrumPtr& input, const double drift_lower, const double drift_upper) - { - OPENMS_PRECONDITION(drift_upper > 0, "Cannot filter by drift time if upper value is less or equal to zero"); - //OPENMS_PRECONDITION(input->getDriftTimeArray() != nullptr, "Cannot filter by drift time if no drift time is available."); - - if (input->getDriftTimeArray() == nullptr) + else // (spectra_addition_method_ == SpectrumAdditionMethod::RESAMPLE) { - std::cerr << "Warning: Cannot filter by drift time if no drift time is available.\n"; - return input; + std::vector spectrum_out; + //added_spec = SpectrumAddition::addUpSpectra(all_spectra, spacing_for_spectra_resampling_, true); + spectrum_out.push_back(SpectrumAddition::addUpSpectra(all_spectra, im_range, spacing_for_spectra_resampling_, true)); + return spectrum_out; } - - OpenSwath::SpectrumPtr output(new OpenSwath::Spectrum); - - OpenSwath::BinaryDataArrayPtr mz_arr = input->getMZArray(); - OpenSwath::BinaryDataArrayPtr int_arr = input->getIntensityArray(); - OpenSwath::BinaryDataArrayPtr im_arr = input->getDriftTimeArray(); - - std::vector::const_iterator mz_it = mz_arr->data.begin(); - std::vector::const_iterator int_it = int_arr->data.begin(); - std::vector::const_iterator im_it = im_arr->data.begin(); - std::vector::const_iterator mz_end = mz_arr->data.end(); - - OpenSwath::BinaryDataArrayPtr mz_arr_out(new OpenSwath::BinaryDataArray); - OpenSwath::BinaryDataArrayPtr intens_arr_out(new OpenSwath::BinaryDataArray); - OpenSwath::BinaryDataArrayPtr im_arr_out(new OpenSwath::BinaryDataArray); - im_arr_out->description = im_arr->description; - - size_t n = mz_arr->data.size(); - im_arr_out->data.reserve(n); - while (mz_it != mz_end) - { - if (*im_it > drift_lower && *im_it < drift_upper) - { - mz_arr_out->data.push_back( *mz_it ); - intens_arr_out->data.push_back( *int_it ); - im_arr_out->data.push_back( *im_it ); - } - ++mz_it; - ++int_it; - ++im_it; - } - output->setMZArray(mz_arr_out); - output->setIntensityArray(intens_arr_out); - output->getDataArrays().push_back(im_arr_out); - return output; } - - OpenSwath::SpectrumPtr OpenSwathScoring::getAddedSpectra_(const OpenSwath::SpectrumAccessPtr& swath_map, - double RT, int nr_spectra_to_add, const double drift_lower, const double drift_upper) + SpectrumSequence OpenSwathScoring::fetchSpectrumSwath(std::vector swath_maps, double RT, int nr_spectra_to_add, const RangeMobility& im_range) { - std::vector indices = swath_map->getSpectraByRT(RT, 0.0); - OpenSwath::SpectrumPtr added_spec(new OpenSwath::Spectrum); - added_spec->getDataArrays().push_back( OpenSwath::BinaryDataArrayPtr(new OpenSwath::BinaryDataArray) ); - added_spec->getDataArrays().back()->description = "Ion Mobility"; + OPENMS_PRECONDITION(nr_spectra_to_add >= 1, "nr_spectra_to_add must be at least 1.") + OPENMS_PRECONDITION(!swath_maps.empty(), "swath_maps vector cannot be empty") - if (indices.empty() ) - { - return added_spec; - } - int closest_idx = boost::numeric_cast(indices[0]); - if (indices[0] != 0 && - std::fabs(swath_map->getSpectrumMetaById(boost::numeric_cast(indices[0]) - 1).RT - RT) < - std::fabs(swath_map->getSpectrumMetaById(boost::numeric_cast(indices[0])).RT - RT)) - { - closest_idx--; - } - - if (nr_spectra_to_add == 1) + // This is not SONAR data + if (swath_maps.size() == 1) { - added_spec = swath_map->getSpectrumById(closest_idx); - if (drift_upper > 0) - { - added_spec = filterByDrift(added_spec, drift_lower, drift_upper); - } + return fetchSpectrumSwath(swath_maps[0].sptr, RT, nr_spectra_to_add, im_range); } else { - std::vector all_spectra; - // always add the spectrum 0, then add those right and left - all_spectra.push_back(swath_map->getSpectrumById(closest_idx)); - for (int i = 1; i <= nr_spectra_to_add / 2; i++) // cast to int is intended! + // data is not IM enhanced + if (!im_range.isEmpty()) { - if (closest_idx - i >= 0) + // multiple SWATH maps for a single precursor -> this is SONAR data, in all cases only return a single spectrum + SpectrumSequence all_spectra; + + if (spectra_addition_method_ == SpectrumAdditionMethod::ADDITION) { - all_spectra.push_back(swath_map->getSpectrumById(closest_idx - i)); + for (size_t i = 0; i < swath_maps.size(); ++i) + { + SpectrumSequence spectrumSequence = swath_maps[i].sptr->getMultipleSpectra(RT, nr_spectra_to_add, im_range.getMin(), im_range.getMax()); + } } - if (closest_idx + i < (int)swath_map->getNrSpectra()) + else // (spectra_addition_method_ == SpectrumAdditionMethod::RESAMPLE) { - all_spectra.push_back(swath_map->getSpectrumById(closest_idx + i)); + for (size_t i = 0; i < swath_maps.size(); ++i) + { + SpectrumSequence spectrumSequence = swath_maps[i].sptr->getMultipleSpectra(RT, nr_spectra_to_add, im_range.getMin(), im_range.getMax()); + all_spectra.push_back(SpectrumAddition::addUpSpectra(spectrumSequence, spacing_for_spectra_resampling_, true)); + } } + return { SpectrumAddition::addUpSpectra(all_spectra, spacing_for_spectra_resampling_, true) }; } - - // Filter all spectra by drift time before further processing - if (drift_upper > 0) + else // im_range.isEmpty() { - for (auto& s: all_spectra) s = filterByDrift(s, drift_lower, drift_upper); - } + // multiple SWATH maps for a single precursor -> this is SONAR data, in all cases only return a single spectrum + SpectrumSequence all_spectra; - // add up all spectra - if (spectra_addition_method_ == "simple") - { - // Ensure that we have the same number of data arrays as in the input spectrum - if (!all_spectra.empty() && all_spectra[0]->getDataArrays().size() > 2) + if (spectra_addition_method_ == SpectrumAdditionMethod::ADDITION) { - for (Size k = 2; k < all_spectra[0]->getDataArrays().size(); k++) + for (size_t i = 0; i < swath_maps.size(); ++i) { - OpenSwath::BinaryDataArrayPtr tmp (new OpenSwath::BinaryDataArray()); - tmp->description = all_spectra[0]->getDataArrays()[k]->description; - added_spec->getDataArrays().push_back(tmp); + SpectrumSequence spectrumSequence = swath_maps[i].sptr->getMultipleSpectra(RT, nr_spectra_to_add); + all_spectra.push_back(SpectrumAddition::concatenateSpectra(spectrumSequence)); } } - - // Simply add up data and sort in the end - for (const auto& s : all_spectra) + else // (spectra_addition_method_ == SpectrumAdditionMethod::RESAMPLE) { - for (Size k = 0; k < s->getDataArrays().size(); k++) + for (size_t i = 0; i < swath_maps.size(); ++i) { - auto& v1 = added_spec->getDataArrays()[k]->data; - auto& v2 = s->getDataArrays()[k]->data; - - v1.reserve( v1.size() + v2.size() ); - v1.insert( v1.end(), v2.begin(), v2.end() ); + SpectrumSequence spectrumSequence = swath_maps[i].sptr->getMultipleSpectra(RT, nr_spectra_to_add); + all_spectra.push_back(SpectrumAddition::addUpSpectra(spectrumSequence, spacing_for_spectra_resampling_, true)); } } - sortSpectrumByMZ(*added_spec); - } - else - { - added_spec = SpectrumAddition::addUpSpectra(all_spectra, spacing_for_spectra_resampling_, true); + return { SpectrumAddition::addUpSpectra(all_spectra, spacing_for_spectra_resampling_, true) }; } } - - OPENMS_POSTCONDITION( std::adjacent_find(added_spec->getMZArray()->data.begin(), - added_spec->getMZArray()->data.end(), std::greater()) == added_spec->getMZArray()->data.end(), - "Postcondition violated: m/z vector needs to be sorted!" ) - - return added_spec; } - } - diff --git a/src/openms/source/ANALYSIS/OPENSWATH/OpenSwathWorkflow.cpp b/src/openms/source/ANALYSIS/OPENSWATH/OpenSwathWorkflow.cpp index d978ab3f6e4..7d13dc7b916 100644 --- a/src/openms/source/ANALYSIS/OPENSWATH/OpenSwathWorkflow.cpp +++ b/src/openms/source/ANALYSIS/OPENSWATH/OpenSwathWorkflow.cpp @@ -230,7 +230,7 @@ namespace OpenMS "There are less than 2 iRT normalization peptides, not enough for an RT correction."); } - // 7. Select the "correct" peaks for m/z correction (e.g. remove those not + // 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) @@ -252,9 +252,11 @@ namespace OpenMS } } - // 8. Correct m/z deviations using SwathMapMassCorrection + // 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); @@ -501,10 +503,6 @@ namespace OpenMS // (i) Obtain precursor chromatograms (MS1) if precursor extraction is enabled ChromExtractParams ms1_cp(cp_ms1); - if (!use_ms1_ion_mobility_) - { - ms1_cp.im_extraction_window = -1; - } if (ms1_only && !use_ms1_traces_) { diff --git a/src/openms/source/ANALYSIS/OPENSWATH/SONARScoring.cpp b/src/openms/source/ANALYSIS/OPENSWATH/SONARScoring.cpp index 903f668ee31..d85982f807a 100644 --- a/src/openms/source/ANALYSIS/OPENSWATH/SONARScoring.cpp +++ b/src/openms/source/ANALYSIS/OPENSWATH/SONARScoring.cpp @@ -187,20 +187,11 @@ namespace OpenMS OpenSwath::SpectrumPtr spectrum_ = swath_map->getSpectrumById(closest_idx); // integrate intensity within that scan - double left = transitions[k].getProductMZ(); - double right = transitions[k].getProductMZ(); - if (dia_extraction_ppm_) - { - left -= left * dia_extract_window_ / 2e6; - right += right * dia_extract_window_ / 2e6; - } - else - { - left -= dia_extract_window_ / 2.0; - right += dia_extract_window_ / 2.0; - } - double mz, intensity; - DIAHelpers::integrateWindow(spectrum_, left, right, mz, intensity, dia_centroided_); + RangeMZ mz_range = DIAHelpers::createMZRangePPM(transitions[k].getProductMZ(), dia_extract_window_, dia_extraction_ppm_); + + double mz, intensity, im; // create im even though not used + RangeMobility im_range; + DIAHelpers::integrateWindow(spectrum_, mz, im, intensity, mz_range, im_range, dia_centroided_); sonar_profile.push_back(intensity); sonar_mz_profile.push_back(mz); diff --git a/src/openms/source/ANALYSIS/OPENSWATH/SpectrumAddition.cpp b/src/openms/source/ANALYSIS/OPENSWATH/SpectrumAddition.cpp index 31fd0dcf77d..1ea4b8947b0 100644 --- a/src/openms/source/ANALYSIS/OPENSWATH/SpectrumAddition.cpp +++ b/src/openms/source/ANALYSIS/OPENSWATH/SpectrumAddition.cpp @@ -11,14 +11,48 @@ #include #include #include +#include // std::iota namespace OpenMS { - OpenSwath::SpectrumPtr SpectrumAddition::addUpSpectra(const std::vector& all_spectra, - double sampling_rate, bool filter_zeros) + void SpectrumAddition::sortSpectrumByMZ(OpenSwath::Spectrum& spec) + { + // Based off of https://stackoverflow.com/questions/1577475/c-sorting-and-keeping-track-of-indexes/43922758#43922758 + + //initialize + //std::vector > sorted_indices; + std::vector sorted_indices(spec.getMZArray()->data.size()); + std::iota(sorted_indices.begin(), sorted_indices.end(), 0); + + // sort indexes based on comparing values in v + // using std::stable_sort instead of std::sort + // to avoid unnecessary index re-orderings + // when v contains elements of equal values + // get the ordering of the indices where mz is sorted + OpenSwath::BinaryDataArrayPtr mzArr = spec.getMZArray(); + std::stable_sort(sorted_indices.begin(), sorted_indices.end(), + [mzArr](size_t i1, size_t i2) {return mzArr->data[i1] < mzArr->data[i2];}); + + // apply sorting across all arrays + for (auto& da : spec.getDataArrays() ) + { + if (da->data.empty()) continue; + for (Size i = 0; i < sorted_indices.size(); ++i) + { + auto j = sorted_indices[i]; + while (jdata[i], da->data[j]); + } + } + + OPENMS_POSTCONDITION( std::adjacent_find(spec.getMZArray()->data.begin(), + spec.getMZArray()->data.end(), std::greater()) == spec.getMZArray()->data.end(), + "Postcondition violated: m/z vector needs to be sorted!" ) + } + + OpenSwath::SpectrumPtr SpectrumAddition::addUpSpectra(const SpectrumSequence& all_spectra, double sampling_rate, bool filter_zeros) { - OPENMS_PRECONDITION(all_spectra.empty() || all_spectra[0]->getDataArrays().size() == 2, "Can only resample spectra with 2 data dimensions (no ion mobility spectra)") if (all_spectra.size() == 1) return all_spectra[0]; if (all_spectra.empty()) @@ -100,6 +134,11 @@ namespace OpenMS if (!filter_zeros) { + OPENMS_POSTCONDITION(std::adjacent_find(resampled_peak_container->getMZArray()->data.begin(), + resampled_peak_container->getMZArray()->data.end(), std::greater()) == resampled_peak_container->getMZArray()->data.end(), + "Postcondition violated: m/z vector needs to be sorted!" ) + + return resampled_peak_container; } else @@ -113,11 +152,77 @@ namespace OpenMS master_spectrum_filtered->getMZArray()->data.push_back(resampled_peak_container->getMZArray()->data[i]); } } + + OPENMS_POSTCONDITION( std::adjacent_find(master_spectrum_filtered->getMZArray()->data.begin(), + master_spectrum_filtered->getMZArray()->data.end(), std::greater()) == master_spectrum_filtered->getMZArray()->data.end(), + "Postcondition violated: m/z vector needs to be sorted!" ) + + return master_spectrum_filtered; } } - OpenMS::MSSpectrum SpectrumAddition::addUpSpectra(const std::vector& all_spectra, double sampling_rate, bool filter_zeros) + + OpenSwath::SpectrumPtr SpectrumAddition::addUpSpectra(const SpectrumSequence& all_spectra, + const RangeMobility& im_range, + double sampling_rate, + bool filter_zeros) + { + OpenSwath::SpectrumPtr added_spec(new OpenSwath::Spectrum); + + // If no spectra found return + if (all_spectra.empty()) + { + return added_spec; + } + + if (im_range.isEmpty()) + { + return addUpSpectra(all_spectra, sampling_rate, filter_zeros); + } + // since resampling is not supported on 3D data first filter by drift time (if possible) and then add + // (!im_range.isEmpty()) + SpectrumSequence filteredSpectra; + for (auto spec: all_spectra) + { + filteredSpectra.push_back(OpenSwath::ISpectrumAccess::filterByDrift(spec, im_range.getMin(), im_range.getMax())); + } + return addUpSpectra(filteredSpectra, sampling_rate, filter_zeros); + } + + OpenSwath::SpectrumPtr SpectrumAddition::concatenateSpectra(const SpectrumSequence& all_spectra) + + { + OpenSwath::SpectrumPtr added_spec(new OpenSwath::Spectrum); + // Ensure that we have the same number of data arrays as in the input spectrum + // copying the extra data arrays descriptions onto the added spectra + if (!all_spectra.empty() && all_spectra[0]->getDataArrays().size() > 2) + { + for (Size k = 2; k < all_spectra[0]->getDataArrays().size(); k++) + { + OpenSwath::BinaryDataArrayPtr tmp (new OpenSwath::BinaryDataArray()); + tmp->description = all_spectra[0]->getDataArrays()[k]->description; + added_spec->getDataArrays().push_back(tmp); + } + } + + // Simply concatenate all spectra together and sort in the end + for (const auto& s : all_spectra) + { + for (Size k = 0; k < s->getDataArrays().size(); k++) + { + auto& v1 = added_spec->getDataArrays()[k]->data; + auto& v2 = s->getDataArrays()[k]->data; + + v1.reserve( v1.size() + v2.size() ); + v1.insert( v1.end(), v2.begin(), v2.end() ); + } + } + sortSpectrumByMZ(*added_spec); + return added_spec; + } + + OpenMS::MSSpectrum SpectrumAddition::addUpSpectra(const std::vector& all_spectra, double sampling_rate, bool filter_zeros) { OPENMS_PRECONDITION(all_spectra.empty() || all_spectra[0].getFloatDataArrays().empty(), "Can only resample spectra with 2 data dimensions (no ion mobility spectra)") @@ -205,6 +310,4 @@ namespace OpenMS return master_spectrum_filtered; } } - } - diff --git a/src/openms/source/ANALYSIS/OPENSWATH/SwathMapMassCorrection.cpp b/src/openms/source/ANALYSIS/OPENSWATH/SwathMapMassCorrection.cpp index 915478b0028..ba8a6984ed7 100644 --- a/src/openms/source/ANALYSIS/OPENSWATH/SwathMapMassCorrection.cpp +++ b/src/openms/source/ANALYSIS/OPENSWATH/SwathMapMassCorrection.cpp @@ -212,27 +212,35 @@ namespace OpenMS #pragma omp critical #endif { + RangeMobility im_range; if (ms1_im_) { - sp_ms1 = OpenSwathScoring().fetchSpectrumSwath(ms1_maps, bestRT, 1, 0, 0); + std::vector fetchSpectrumArr = OpenSwathScoring().fetchSpectrumSwath(ms1_maps, bestRT, 1, im_range); + sp_ms1 = (!fetchSpectrumArr.empty()) ? fetchSpectrumArr[0] : *new(OpenSwath::SpectrumPtr); } else { - sp_ms2 = OpenSwathScoring().fetchSpectrumSwath(used_maps, bestRT, 1, 0, 0); + std::vector fetchSpectrumArr = OpenSwathScoring().fetchSpectrumSwath(used_maps, bestRT, 1, im_range); + sp_ms2 = (!fetchSpectrumArr.empty()) ? fetchSpectrumArr[0] : *new(OpenSwath::SpectrumPtr); } } for (const auto& tr : transition_group->getTransitions()) { if (ms1_im_) {continue;} - double intensity(0), im(0), left(tr.product_mz), right(tr.product_mz); + double intensity(0), im(0), mz(0); + RangeMZ mz_range = DIAHelpers::createMZRangePPM(tr.product_mz, mz_extr_window, ppm); // get drift time upper/lower offset (this assumes that all chromatograms // are derived from the same precursor with the same drift time) auto pepref = tr.getPeptideRef(); double drift_target = pep_im_map[pepref]; - double drift_left(drift_target), drift_right(drift_target); - DIAHelpers::adjustExtractionWindow(drift_right, drift_left, im_extraction_win, false); + RangeMobility im_range; + if (im_extraction_win != -1 ) // im_extraction_win is set + { + im_range.setMax(drift_target); + im_range.minSpanIfSingular(im_extraction_win); + } // Check that the spectrum really has a drift time array if (sp_ms2->getDriftTimeArray() == nullptr) @@ -245,8 +253,7 @@ namespace OpenMS continue; } - DIAHelpers::adjustExtractionWindow(right, left, mz_extr_window, ppm); - DIAHelpers::integrateDriftSpectrum(sp_ms2, left, right, im, intensity, drift_left, drift_right); + DIAHelpers::integrateWindow(sp_ms2, mz, im, intensity, mz_range, im_range); // skip empty windows if (im <= 0) @@ -274,14 +281,17 @@ namespace OpenMS if (!transition_group->getTransitions().empty() && ms1_im_) { const auto& tr = transition_group->getTransitions()[0]; - double intensity(0), im(0), left(tr.precursor_mz), right(tr.precursor_mz); + double intensity(0), im(0), mz(0); + RangeMZ mz_range = DIAHelpers::createMZRangePPM(tr.precursor_mz, mz_extr_window, ppm); // get drift time upper/lower offset (this assumes that all chromatograms // are derived from the same precursor with the same drift time) auto pepref = tr.getPeptideRef(); double drift_target = pep_im_map[pepref]; - double drift_left(drift_target), drift_right(drift_target); - DIAHelpers::adjustExtractionWindow(drift_right, drift_left, im_extraction_win, false); + + // do not need to check for IM because we are correcting IM + RangeMobility im_range(drift_target); + im_range.minSpanIfSingular(im_extraction_win); // Check that the spectrum really has a drift time array if (sp_ms1->getDriftTimeArray() == nullptr) @@ -294,8 +304,7 @@ namespace OpenMS continue; } - DIAHelpers::adjustExtractionWindow(right, left, mz_extr_window, ppm); - DIAHelpers::integrateDriftSpectrum(sp_ms1, left, right, im, intensity, drift_left, drift_right); + DIAHelpers::integrateWindow(sp_ms1, mz, im, intensity, mz_range, im_range); // skip empty windows if (im <= 0) @@ -358,6 +367,7 @@ namespace OpenMS bool ppm = mz_extraction_window_ppm_; double mz_extr_window = mz_extraction_window_; std::string corr_type = mz_correction_function_; + double im_extraction = im_extraction_window_; OPENMS_LOG_DEBUG << "SwathMapMassCorrection::correctMZ with type " << corr_type << " and window " << mz_extr_window << " in ppm " << ppm << std::endl; @@ -422,17 +432,26 @@ namespace OpenMS continue; } + // if ion mobility extraction window is set than extract with ion mobility + RangeMobility im_range; + if (im_extraction != -1) // ion mobility extraction is set + { + im_range.setMax(drift_target); + im_range.minSpanIfSingular(im_extraction); + } + // Get the spectrum for this RT and extract raw data points for all the // calibrating transitions (fragment m/z values) from the spectrum - OpenSwath::SpectrumPtr sp = OpenSwathScoring().fetchSpectrumSwath(used_maps, bestRT, 1, 0, 0); + std::vector spArr = OpenSwathScoring().fetchSpectrumSwath(used_maps, bestRT, 1, im_range); + OpenSwath::SpectrumPtr sp = (!spArr.empty()) ? spArr[0] : *new(OpenSwath::SpectrumPtr); for (const auto& tr : transition_group->getTransitions()) { - double mz, intensity, left(tr.product_mz), right(tr.product_mz); + double mz, intensity, im; + RangeMZ mz_range = DIAHelpers::createMZRangePPM(tr.product_mz, mz_extr_window, ppm); bool centroided = false; // integrate spectrum at the position of the theoretical mass - DIAHelpers::adjustExtractionWindow(right, left, mz_extr_window, ppm); - DIAHelpers::integrateWindow(sp, left, right, mz, intensity, centroided); + DIAHelpers::integrateWindow(sp, mz, im, intensity, mz_range, im_range, centroided); // Correct using the irt_im // skip empty windows if (mz == -1) diff --git a/src/openswathalgo/include/OpenMS/OPENSWATHALGO/DATAACCESS/DataStructures.h b/src/openswathalgo/include/OpenMS/OPENSWATHALGO/DATAACCESS/DataStructures.h index 36e3c3729b2..d185953bb9e 100644 --- a/src/openswathalgo/include/OpenMS/OPENSWATHALGO/DATAACCESS/DataStructures.h +++ b/src/openswathalgo/include/OpenMS/OPENSWATHALGO/DATAACCESS/DataStructures.h @@ -231,6 +231,12 @@ namespace OpenSwath binaryDataArrayPtrs[1] = data; } + void setDriftTimeArray(BinaryDataArrayPtr data) + { + data->description = "Ion Mobility"; + binaryDataArrayPtrs.push_back(data); + } + /// get drift time array (may be null) BinaryDataArrayPtr getDriftTimeArray() const { diff --git a/src/openswathalgo/include/OpenMS/OPENSWATHALGO/DATAACCESS/ISpectrumAccess.h b/src/openswathalgo/include/OpenMS/OPENSWATHALGO/DATAACCESS/ISpectrumAccess.h index b475b8add79..baa306ad8a5 100644 --- a/src/openswathalgo/include/OpenMS/OPENSWATHALGO/DATAACCESS/ISpectrumAccess.h +++ b/src/openswathalgo/include/OpenMS/OPENSWATHALGO/DATAACCESS/ISpectrumAccess.h @@ -15,9 +15,14 @@ #include #include +namespace OpenMS +{ + using SpectrumSequence = std::vector; ///< a vector of spectrum pointers that DIA scores can operate on, allows for clever integration of only the target regions +} namespace OpenSwath { + using SpectrumSequence = OpenMS::SpectrumSequence; /** @brief The interface of a mass spectrometry experiment. */ @@ -43,6 +48,10 @@ namespace OpenSwath /// Return a pointer to a spectrum at the given id virtual SpectrumPtr getSpectrumById(int id) = 0; + + /// Return pointer to a spectrum at the given id, the spectrum will be filtered by drift time + SpectrumPtr getSpectrumById(int id, double drift_start, double drift_end ); + /// Return a vector of ids of spectra that are within RT +/- deltaRT virtual std::vector getSpectraByRT(double RT, double deltaRT) const = 0; /// Returns the number of spectra available @@ -56,7 +65,68 @@ namespace OpenSwath virtual std::size_t getNrChromatograms() const = 0; /// Returns the native id of the chromatogram at the given id virtual std::string getChromatogramNativeID(int id) const = 0; - }; + + /* @brief Fetches a spectrumSequence (multiple spectra pointers) closest to the given RT + * @p RT = target RT + * @p nr_spectra_to_fetch = # spectra around target RT to fetch (length of the spectrum sequence) + */ + SpectrumSequence getMultipleSpectra(double RT, int nr_spectra_to_fetch); + + /* @brief Fetches a spectrumSequence (multiple spectra pointers) closest to the given RT. Filters all spectra by specified @p drift_start and @p drift_end + * @p RT = target RT + * @p nr_spectra_to_fetch = # spectra around target RT to fetch (length of the spectrum sequence) + */ + SpectrumSequence getMultipleSpectra(double RT, int nr_spectra_to_fetch, double drift_start, double drift_end); + + /// filters a spectrum by drift time, spectrum pointer returned is a copy + static SpectrumPtr filterByDrift(const SpectrumPtr& input, double drift_start, double drift_end) + { + // NOTE: this function is very inefficient because filtering unsorted array + //OPENMS_PRECONDITION(drift_start <= 0, "Cannot filter by drift time if drift_start is not set"); + //OPENMS_PRECONDITION(drift_end - drift_start < 0, "Cannot filter by drift time if range is empty"); + //OPENMS_PRECONDITION(input->getDriftTimeArray() != nullptr, "Cannot filter by drift time if no drift time is available."); + + //if (input->getDriftTimeArray() == nullptr) + //{ + //throw Exception::NullPointer(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION); + //} + + OpenSwath::SpectrumPtr output(new OpenSwath::Spectrum); + + OpenSwath::BinaryDataArrayPtr mz_arr = input->getMZArray(); + OpenSwath::BinaryDataArrayPtr int_arr = input->getIntensityArray(); + OpenSwath::BinaryDataArrayPtr im_arr = input->getDriftTimeArray(); + + auto mz_it = mz_arr->data.cbegin(); + auto int_it = int_arr->data.cbegin(); + auto im_it = im_arr->data.cbegin(); + auto mz_end = mz_arr->data.cend(); + + OpenSwath::BinaryDataArrayPtr mz_arr_out(new OpenSwath::BinaryDataArray); + OpenSwath::BinaryDataArrayPtr intens_arr_out(new OpenSwath::BinaryDataArray); + OpenSwath::BinaryDataArrayPtr im_arr_out(new OpenSwath::BinaryDataArray); + im_arr_out->description = im_arr->description; + + while (mz_it != mz_end) + { + if ( (drift_start <= *im_it) & (drift_end >= *im_it) ) + { + mz_arr_out->data.push_back( *mz_it ); + intens_arr_out->data.push_back( *int_it ); + im_arr_out->data.push_back( *im_it ); + } + ++mz_it; + ++int_it; + ++im_it; + } + output->setMZArray(mz_arr_out); + output->setIntensityArray(intens_arr_out); + output->getDataArrays().push_back(im_arr_out); + return output; + } + + + }; typedef boost::shared_ptr SpectrumAccessPtr; } diff --git a/src/openswathalgo/source/OPENSWATHALGO/DATAACCESS/ISpectrumAccess.cpp b/src/openswathalgo/source/OPENSWATHALGO/DATAACCESS/ISpectrumAccess.cpp index 4f628cdeafe..f018da79648 100644 --- a/src/openswathalgo/source/OPENSWATHALGO/DATAACCESS/ISpectrumAccess.cpp +++ b/src/openswathalgo/source/OPENSWATHALGO/DATAACCESS/ISpectrumAccess.cpp @@ -8,6 +8,10 @@ #include +#include +#include +#include +#include namespace OpenSwath { @@ -15,4 +19,84 @@ namespace OpenSwath { } + SpectrumSequence ISpectrumAccess::getMultipleSpectra(double RT, int nr_spectra_to_fetch) + { + std::vector indices = getSpectraByRT(RT, 0.0); + SpectrumSequence all_spectra; + + if (indices.empty() ) + { + return all_spectra; + } + int closest_idx = boost::numeric_cast(indices[0]); + if (indices[0] != 0 && + std::fabs(getSpectrumMetaById(boost::numeric_cast(indices[0]) - 1).RT - RT) < + std::fabs(getSpectrumMetaById(boost::numeric_cast(indices[0])).RT - RT)) + { + closest_idx--; + } + + all_spectra.push_back(getSpectrumById(closest_idx)); + + int nrSpectra = (int) getNrSpectra(); + for (int i = 1; i <= nr_spectra_to_fetch / 2; i++) // cast to int is intended! + { + if (closest_idx - i >= 0) + { + all_spectra.push_back(getSpectrumById(closest_idx - i)); + } + if (closest_idx + i < nrSpectra) + { + all_spectra.push_back(getSpectrumById(closest_idx + i)); + } + } + + return all_spectra; + } + + + SpectrumSequence ISpectrumAccess::getMultipleSpectra(double RT, int nr_spectra_to_fetch, double drift_start, double drift_end) + { + std::vector indices = getSpectraByRT(RT, 0.0); + SpectrumSequence all_spectra; + + if (indices.empty() ) + { + return all_spectra; + } + int closest_idx = boost::numeric_cast(indices[0]); + if (indices[0] != 0 && + std::fabs(getSpectrumMetaById(boost::numeric_cast(indices[0]) - 1).RT - RT) < + std::fabs(getSpectrumMetaById(boost::numeric_cast(indices[0])).RT - RT)) + { + closest_idx--; + } + + all_spectra.push_back(getSpectrumById(closest_idx, drift_start, drift_end)); + + int nrSpectra = (int) getNrSpectra(); + for (int i = 1; i <= nr_spectra_to_fetch / 2; i++) // cast to int is intended! + { + if (closest_idx - i >= 0) + { + all_spectra.push_back(getSpectrumById(closest_idx - i, drift_start, drift_end)); + } + if (closest_idx + i < nrSpectra) + { + all_spectra.push_back(getSpectrumById(closest_idx + i, drift_start, drift_end)); + } + } + + return all_spectra; + } + + + SpectrumPtr ISpectrumAccess::getSpectrumById(int id, double drift_start, double drift_end) + { + // first fetch the spectrum + OpenSwath::SpectrumPtr spectrum = getSpectrumById(id); + + // then filter by drift + return ISpectrumAccess::filterByDrift(spectrum, drift_start, drift_end); + } } diff --git a/src/pyOpenMS/addons/DIAScoring.pyx b/src/pyOpenMS/addons/DIAScoring.pyx index 44cb67c9879..b9bc32e19dc 100644 --- a/src/pyOpenMS/addons/DIAScoring.pyx +++ b/src/pyOpenMS/addons/DIAScoring.pyx @@ -1,18 +1,22 @@ - def dia_by_ion_score(self, OSSpectrum spectrum , AASequence sequence , charge , float bseries_score , float yseries_score ): - assert isinstance(spectrum, OSSpectrum), 'arg spectrum wrong type' + def dia_by_ion_score(self, list spectrum , AASequence sequence, charge, RangeMobility im_range, float bseries_score , float yseries_score): + assert isinstance(spectrum, list) and all(isinstance(elemt_rec, OSSpectrum) for elemt_rec in spectrum), 'arg spectrum wrong type' assert isinstance(sequence, AASequence), 'arg sequence wrong type' + assert isinstance(im_range, RangeMobility), 'arg sequence wrong type' assert isinstance(charge, (int, long)), 'arg charge wrong type' assert isinstance(bseries_score, float), 'arg bseries_score wrong type' assert isinstance(yseries_score, float), 'arg yseries_score wrong type' - cdef shared_ptr[_OSSpectrum] input_spectrum = spectrum.inst + + cdef libcpp_vector[shared_ptr[_OSSpectrum]] v1 + cdef OSSpectrum spectrum_rec + for spectrum_rec in spectrum: + v1.push_back(spectrum_rec.inst) cdef double input_bseries_score = (bseries_score) cdef double input_yseries_score = (yseries_score) - self.inst.get().dia_by_ion_score(input_spectrum, (deref(sequence.inst.get())), (charge), input_bseries_score, input_yseries_score) + self.inst.get().dia_by_ion_score(v1, (deref(sequence.inst.get())), (charge), (deref(im_range.inst.get())), input_bseries_score, input_yseries_score) yseries_score = input_yseries_score bseries_score = input_bseries_score return (bseries_score,yseries_score) - diff --git a/src/pyOpenMS/pxds/DIAScoring.pxd b/src/pyOpenMS/pxds/DIAScoring.pxd index ede49031962..d81f1555809 100644 --- a/src/pyOpenMS/pxds/DIAScoring.pxd +++ b/src/pyOpenMS/pxds/DIAScoring.pxd @@ -6,6 +6,8 @@ from LightTargetedExperiment cimport * from AASequence cimport * from DefaultParamHandler cimport * from EmpiricalFormula cimport * +from libcpp.vector cimport vector as libcpp_vector +from RangeManager cimport * cdef extern from "" namespace "OpenMS": @@ -17,18 +19,17 @@ cdef extern from "" namespace "OpenMS": # private DIAScoring(DIAScoring) except + nogil # wrap-ignore - bool dia_ms1_massdiff_score(double precursor_mz, OSSpectrumPtr spectrum, - double& ppm_score) except + nogil + bool dia_ms1_massdiff_score(double precursor_mz, libcpp_vector[shared_ptr[OSSpectrum] ] spectrum, + RangeMobility& im_range, double& ppm_score) except + nogil - void dia_ms1_isotope_scores_averagine(double precursor_mz, OSSpectrumPtr spectrum, - double& isotope_corr, double& isotope_overlap, int charge_state) except + nogil + void dia_ms1_isotope_scores_averagine(double precursor_mz, libcpp_vector[shared_ptr[OSSpectrum] ] spectrum, int charge_state, RangeMobility& im_range, + double& isotope_corr, double& isotope_overlap) except + nogil - void dia_ms1_isotope_scores(double precursor_mz, OSSpectrumPtr spectrum, - double& isotope_corr, double& isotope_overlap, EmpiricalFormula& sum_formula) except + nogil + void dia_ms1_isotope_scores(double precursor_mz, libcpp_vector[shared_ptr[OSSpectrum] ] spectrum, RangeMobility& im_range, + double& isotope_corr, double& isotope_overlap, EmpiricalFormula& sum_formula) except + nogil # TODO automatically wrap - void dia_by_ion_score(OSSpectrumPtr spectrum, AASequence sequence, int charge, double & bseries_score, double & yseries_score) except + nogil # wrap-return:return(bseries_score,yseries_score) wrap-ignore + void dia_by_ion_score( libcpp_vector[shared_ptr[OSSpectrum] ] spectrum, AASequence sequence, int charge, RangeMobility& im_range, double & bseries_score, double & yseries_score) except + nogil # wrap-return:return(bseries_score,yseries_score) wrap-ignore # Dotproduct / Manhatten score with theoretical spectrum - void score_with_isotopes(OSSpectrumPtr spectrum, libcpp_vector[LightTransition] transitions, - double& dotprod, double& manhattan) except + nogil - + void score_with_isotopes( libcpp_vector[shared_ptr[OSSpectrum] ] spectrum, libcpp_vector[LightTransition] transitions, RangeMobility& im_range, + double& dotprod, double& manhattan) except + nogil diff --git a/src/pyOpenMS/pxds/OpenSwathScoring.pxd b/src/pyOpenMS/pxds/OpenSwathScoring.pxd index da8a0880c94..e6051ab349d 100644 --- a/src/pyOpenMS/pxds/OpenSwathScoring.pxd +++ b/src/pyOpenMS/pxds/OpenSwathScoring.pxd @@ -16,7 +16,8 @@ cdef extern from "" namespace "Ope double spacing_for_spectra_resampling, double drift_extra, OpenSwath_Scores_Usage su, - libcpp_string spectrum_addition_method) except + nogil + libcpp_string spectrum_addition_method, + bool use_ms1_ion_mobility) except + nogil # wrap-doc: # Initialize the scoring object\n # Sets the parameters for the scoring diff --git a/src/pyOpenMS/pxds/RangeManager.pxd b/src/pyOpenMS/pxds/RangeManager.pxd index af8ebeef1b2..9901c434770 100644 --- a/src/pyOpenMS/pxds/RangeManager.pxd +++ b/src/pyOpenMS/pxds/RangeManager.pxd @@ -31,10 +31,40 @@ cdef extern from "" namespace "OpenMS": cdef cppclass RangeManagerRtInt "OpenMS::RangeManager": # wrap-ignore # no-pxd-import - RangeManagerRtInt() except + nogil - RangeManagerRtInt(RangeManagerRtInt &) except + nogil - double getMinRT() except + nogil # wrap-doc:Returns the minimum RT - double getMaxRT() except + nogil # wrap-doc:Returns the maximum RT - double getMinIntensity() except + nogil # wrap-doc:Returns the minimum intensity - double getMaxIntensity() except + nogil # wrap-doc:Returns the maximum intensity - void clearRanges() except + nogil # wrap-doc:Resets all range dimensions as empty + RangeManagerRtInt() except + nogil + RangeManagerRtInt(RangeManagerRtInt &) except + nogil + double getMinRT() except + nogil # wrap-doc:Returns the minimum RT + double getMaxRT() except + nogil # wrap-doc:Returns the maximum RT + double getMinIntensity() except + nogil # wrap-doc:Returns the minimum intensity + double getMaxIntensity() except + nogil # wrap-doc:Returns the maximum intensity + void clearRanges() except + nogil # wrap-doc:Resets all range dimensions as empty + + cdef cppclass RangeBase: + # wrap-ignore + # no-pxd-import + + RangeBase() except + nogil + RangeBase(RangeBase &) except + nogil + + cdef cppclass RangeMZ (RangeBase) : + # wrap-inherits: + # RangeBase + + RangeMZ() except + nogil + RangeMZ(RangeMZ &) except + nogil #compiler + + + cdef cppclass RangeIntensity (RangeBase) : + # wrap-inherits: + # RangeBase + + RangeIntensity() except + nogil + RangeIntensity(RangeIntensity &) except + nogil + + + cdef cppclass RangeMobility (RangeBase) : + # wrap-inherits: + # RangeBase + + RangeMobility() except + nogil + RangeMobility(RangeMobility &) except + nogil diff --git a/src/pyOpenMS/tests/unittests/test_DIAScoring.py b/src/pyOpenMS/tests/unittests/test_DIAScoring.py index 5ef91b3fb0d..dafcb6c0dca 100644 --- a/src/pyOpenMS/tests/unittests/test_DIAScoring.py +++ b/src/pyOpenMS/tests/unittests/test_DIAScoring.py @@ -33,6 +33,9 @@ def test_spectrum(self): spectrum.setMZArray(mz) spectrum.setIntensityArray(intensity) + + spectrumList = [ spectrum ] + diascoring = pyopenms.DIAScoring() # diascoring.set_dia_parameters(0.05, False, 30, 50, 4, 4) // here we use a large enough window so that none of our peaks falls out p_dia = diascoring.getDefaults(); @@ -50,7 +53,8 @@ def test_spectrum(self): bseries_score = 0.0 yseries_score = 0.0 charge = 1 - bseries_score, yseries_score = diascoring.dia_by_ion_score(spectrum, a, charge, bseries_score, yseries_score) + im_range = pyopenms.RangeMobility() + bseries_score, yseries_score = diascoring.dia_by_ion_score([spectrum], a, charge, im_range, bseries_score, yseries_score) self.assertAlmostEqual(bseries_score, 2.0) self.assertAlmostEqual(yseries_score, 2.0) @@ -59,7 +63,8 @@ def test_spectrum(self): a.setModification(1, b"Phospho" ) #; // modify the Y bseries_score = 0 yseries_score = 0 - bseries_score, yseries_score = diascoring.dia_by_ion_score(spectrum, a, 1, bseries_score, yseries_score) + im_range = pyopenms.RangeMobility() + bseries_score, yseries_score = diascoring.dia_by_ion_score([spectrum], a, 1, im_range, bseries_score, yseries_score) self.assertAlmostEqual (bseries_score, 1.0) self.assertAlmostEqual (yseries_score, 3.0) diff --git a/src/tests/class_tests/openms/source/DIAHelper_test.cpp b/src/tests/class_tests/openms/source/DIAHelper_test.cpp index e5d9e184b73..de98a6f5fd0 100644 --- a/src/tests/class_tests/openms/source/DIAHelper_test.cpp +++ b/src/tests/class_tests/openms/source/DIAHelper_test.cpp @@ -25,76 +25,437 @@ START_TEST(DIAHelper, "$Id$") ///////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////// +OpenSwath::SpectrumPtr imSpec(new OpenSwath::Spectrum()); +{ + OpenSwath::BinaryDataArrayPtr mass(new OpenSwath::BinaryDataArray); + mass->data.push_back(100.); + mass->data.push_back(101.); + mass->data.push_back(102.); + mass->data.push_back(103.); + mass->data.push_back(104.); + mass->data.push_back(105.); + mass->data.push_back(106.); + + OpenSwath::BinaryDataArrayPtr intensity(new OpenSwath::BinaryDataArray); + intensity->data.push_back(1.); + intensity->data.push_back(1.); + intensity->data.push_back(1.); + intensity->data.push_back(1.); + intensity->data.push_back(1.); + intensity->data.push_back(1.); + intensity->data.push_back(1.); + + OpenSwath::BinaryDataArrayPtr ion_mobility(new OpenSwath::BinaryDataArray); + ion_mobility->data.push_back(1.); + ion_mobility->data.push_back(2.); + ion_mobility->data.push_back(3.); + ion_mobility->data.push_back(4.); + ion_mobility->data.push_back(5.); + ion_mobility->data.push_back(6.); + ion_mobility->data.push_back(7.); + ion_mobility->description = "Ion Mobility"; + + imSpec->setMZArray( mass ); + imSpec->setIntensityArray( intensity ); + imSpec->getDataArrays().push_back( ion_mobility ); +} + +OpenSwath::SpectrumPtr spec(new OpenSwath::Spectrum()); +{ + OpenSwath::BinaryDataArrayPtr mass(new OpenSwath::BinaryDataArray); + mass->data.push_back(100.); + mass->data.push_back(101.); + mass->data.push_back(102.); + mass->data.push_back(103.); + mass->data.push_back(104.); + mass->data.push_back(105.); + mass->data.push_back(106.); + + OpenSwath::BinaryDataArrayPtr intensity(new OpenSwath::BinaryDataArray); + intensity->data.push_back(1.); + intensity->data.push_back(1.); + intensity->data.push_back(1.); + intensity->data.push_back(1.); + intensity->data.push_back(1.); + intensity->data.push_back(1.); + intensity->data.push_back(1.); + + spec->setMZArray( mass ); + spec->setIntensityArray( intensity ); +} + + +START_SECTION(bool integrateWindow(const OpenSwath::SpectrumPtr& spectrum, double & mz, double & im, double & intensity, RangeMZ range_mz, RangeMobility im_range, bool centroided)) +{ + + RangeMobility empty_im_range; // empty im range for testing + + // IM range from 2 to 5 + RangeMobility nonempty_im_range(3.5); + nonempty_im_range.minSpanIfSingular(3); + + //Mz range from 101 to 103 + RangeMZ mz_range(102.); + mz_range.minSpanIfSingular(2.); // not in ppm + // + //mz range from 101 to 109 + RangeMZ mz_range_2(105.); + mz_range_2.minSpanIfSingular(8.); // not in ppm + + { + //Test integration of empty spectrum + OpenSwath::SpectrumPtr emptySpec(new OpenSwath::Spectrum()); + double mz(0), intens(0), im(0); + + DIAHelpers::integrateWindow(emptySpec, mz, im, intens, mz_range, empty_im_range); + + TEST_REAL_SIMILAR(mz, -1); + TEST_REAL_SIMILAR(im, -1); + TEST_REAL_SIMILAR(intens, 0); + } + + { + // Test spectrum without ion mobility while asking for ion mobility filtering, should throw an exception + double mz(0), intens(0), im(0); + + // unfortunately TEST_EXCEPTION_WITH_MESSAGE was failing for me so test without + TEST_EXCEPTION(Exception::MissingInformation, DIAHelpers::integrateWindow(spec, mz, im, intens, mz_range, nonempty_im_range)); + //TEST_EXCEPTION_WITH_MESSAGE(Exception::MissingInformation, DIAHelpers::integrateWindow(spec, mz, im, intens, mz_range, nonempty_im_range), "a"); + // + //TEST_EXCEPTION_WITH_MESSAGE(Exception::MissingInformation, DIAHelpers::integrateWindow(spec, mz, im, intens, mz_range, nonempty_im_range), "Cannot integrate with drift time if no drift time is available"); + } + + { + // Test ion mobility enhanced array with no ion mobility windows, although IM is present it should be ignored + double mz(0), intens(0), im(0); + + DIAHelpers::integrateWindow(imSpec, mz, im, intens, mz_range, empty_im_range); + TEST_REAL_SIMILAR (mz, 101.5); + TEST_REAL_SIMILAR (intens, 2); + TEST_REAL_SIMILAR (im, -1); // since no IM, this value should be -1 + } + + + + { + // Test With Ion Mobility (Condition 1/2) + double mz(0), intens(0), im(0); + + DIAHelpers::integrateWindow(imSpec, mz, im, intens, mz_range_2, nonempty_im_range); + TEST_REAL_SIMILAR (im, 3.5); + TEST_REAL_SIMILAR (intens, 4); + } + + { + // Test with Ion Mobility (Condition 2/2) + double mz(0), intens(0), im(0); + + DIAHelpers::integrateWindow(imSpec, mz, im, intens, mz_range, nonempty_im_range); + TEST_REAL_SIMILAR (im, 2.5); + TEST_REAL_SIMILAR (intens, 2); + } +} +END_SECTION + + +START_SECTION(bool integrateWindow(const SpectrumSequence& spectra, double & mz, double & im, double & intensity, RangeMZ mz_range, const RangeMobility& im_range, bool centroided)) +{ + + RangeMobility empty_im_range; // empty im range for testing + + // IM range from 2 to 5 + RangeMobility nonempty_im_range(3.5); + nonempty_im_range.minSpanIfSingular(3); + + //Mz range from 101 to 103 + RangeMZ mz_range(102.); // not in ppm + mz_range.minSpanIfSingular(2.); + + + //mz range from 101 to 109 + RangeMZ mz_range_2(105.); + mz_range_2.minSpanIfSingular(8.); // not in ppm + + + { + // Test integration of empty array + std::vector emptySpecArr; + double mz(0), intens(0), im(0); + + DIAHelpers::integrateWindow(emptySpecArr, mz, im, intens, mz_range, empty_im_range); + TEST_REAL_SIMILAR(mz, -1); + TEST_REAL_SIMILAR(im, -1); + TEST_REAL_SIMILAR(intens, 0); + } + + { + //Test integration of empty spectrum + OpenSwath::SpectrumPtr emptySpec(new OpenSwath::Spectrum()); + std::vector specArrEmptySpectrum; + double mz(0), intens(0), im(0); + + specArrEmptySpectrum.push_back(emptySpec); + DIAHelpers::integrateWindow(specArrEmptySpectrum, mz, im, intens, mz_range, empty_im_range); + + TEST_REAL_SIMILAR(mz, -1); + TEST_REAL_SIMILAR(im, -1); + TEST_REAL_SIMILAR(intens, 0); + } + + { + // Test ion mobility enhanced array with no ion mobility windows, although IM is present it should be ignored + std::vector specArr; + double mz(0), intens(0), im(0); + specArr.push_back(imSpec); + + DIAHelpers::integrateWindow(specArr, mz, im, intens, mz_range, empty_im_range); + TEST_REAL_SIMILAR (mz, 101.5); + TEST_REAL_SIMILAR (intens, 2); + TEST_REAL_SIMILAR (im, -1); // since no IM, this value should be -1 + } -START_SECTION(testIntegrateWindows_test) + { + // Test With Ion Mobility (Condition 1/2) + std::vector specArr; + double mz(0), intens(0), im(0); + + specArr.push_back(imSpec); + + DIAHelpers::integrateWindow(specArr, mz, im, intens, mz_range_2, nonempty_im_range); + TEST_REAL_SIMILAR (im, 3.5); + TEST_REAL_SIMILAR (intens, 4); + } + + { + // Test with Ion Mobility (Condition 2/2) + std::vector specArr; + double mz(0), intens(0), im(0); + + specArr.push_back(imSpec); + + DIAHelpers::integrateWindow(specArr, mz, im, intens, mz_range, nonempty_im_range); + TEST_REAL_SIMILAR (im, 2.5); + TEST_REAL_SIMILAR (intens, 2); + } +} +END_SECTION + +START_SECTION(void integrateWindows(const OpenSwath::SpectrumPtr& spectra, const std::vector & windowsCenter, double width, std::vector & integratedWindowsIntensity, std::vector & integratedWindowsMZ, std::vector & integratedWindowsIm, const RangeMobility& im_range, bool remZero)) { - OpenSwath::SpectrumPtr spec(new OpenSwath::Spectrum()); - OpenSwath::BinaryDataArrayPtr mass(new OpenSwath::BinaryDataArray); - - mass->data.push_back(100.); - mass->data.push_back(101.); - mass->data.push_back(102.); - mass->data.push_back(103.); - mass->data.push_back(104.); - mass->data.push_back(105.); - mass->data.push_back(106.); - - OpenSwath::BinaryDataArrayPtr intensity(new OpenSwath::BinaryDataArray); - intensity->data.push_back(1.); - intensity->data.push_back(1.); - intensity->data.push_back(1.); - intensity->data.push_back(1.); - intensity->data.push_back(1.); - intensity->data.push_back(1.); - intensity->data.push_back(1.); - - OpenSwath::BinaryDataArrayPtr ion_mobility(new OpenSwath::BinaryDataArray); - ion_mobility->data.push_back(1.); - ion_mobility->data.push_back(2.); - ion_mobility->data.push_back(3.); - ion_mobility->data.push_back(4.); - ion_mobility->data.push_back(5.); - ion_mobility->data.push_back(6.); - ion_mobility->data.push_back(7.); - ion_mobility->description = "Ion Mobility"; - - spec->setMZArray( mass); - spec->setIntensityArray( intensity); - spec->getDataArrays().push_back( ion_mobility ); - - double mz, intens; - DIAHelpers::integrateWindow(spec, 101., 103., mz, intens); - // std::cout << "mz : " << mz << " int : " << intens << std::endl; - TEST_REAL_SIMILAR (mz, 101.5); - TEST_REAL_SIMILAR (intens, 2) - - std::vector windows, intInt, intMz; - windows.push_back(101.); - windows.push_back(103.); - windows.push_back(105.); - DIAHelpers::integrateWindows(spec, windows, 2, intInt, intMz); - TEST_EQUAL (intInt.size(), intMz.size() ) - TEST_EQUAL (intInt.size(), 3) - TEST_REAL_SIMILAR (intInt[0], 2) - TEST_REAL_SIMILAR (intMz[0], 100.5); - - // std::cout << "print Int" << std::endl; - // std::copy(intInt.begin(), intInt.end(), - // std::ostream_iterator(std::cout, " ")); - // std::cout << std::endl << "print mz" << intMz.size() << std::endl; - // std::cout << intMz[0] << " " << intMz[1] << " " << intMz[2] << std::endl; - // std::copy(intMz.begin(), intMz.end(), - // std::ostream_iterator(std::cout, " ")); - - double im(0.0), im_intens(0.0); - DIAHelpers::integrateDriftSpectrum(spec, 101., 109., im, im_intens, 2, 5); - TEST_REAL_SIMILAR (im, 3.5); - TEST_REAL_SIMILAR (im_intens, 4) - - double im2(0.0), im_intens2(0.0); - DIAHelpers::integrateDriftSpectrum(spec, 101., 103., im2, im_intens2, 2, 5); - TEST_REAL_SIMILAR (im2, 2.5); - TEST_REAL_SIMILAR (im_intens2, 2) + + RangeMobility empty_im_range; // empty im range for testing + RangeMobility nonempty_im_range(3.5); + nonempty_im_range.minSpanIfSingular(3); + + { + // Test empty spectrum (with non empty windows) - remove zeros + OpenSwath::SpectrumPtr emptySpec(new OpenSwath::Spectrum()); + std::vector windows, intInt, intMz, intIm; + + windows.push_back(101.); + windows.push_back(103.); + windows.push_back(105.); + + DIAHelpers::integrateWindows(emptySpec, windows, 2, intInt, intMz, intIm, empty_im_range, true); + TEST_EQUAL (intInt.empty(), true); + TEST_EQUAL (intIm.empty(), true); + TEST_EQUAL (intMz.empty(), true); + } + + { + // Test empty spectrum (with non empty windows) - Don't remove zeros + OpenSwath::SpectrumPtr emptySpec(new OpenSwath::Spectrum()); + std::vector specArr; + std::vector windows, intInt, intMz, intIm; + + windows.push_back(101.); + windows.push_back(103.); + windows.push_back(105.); + + DIAHelpers::integrateWindows(emptySpec, windows, 2, intInt, intMz, intIm, empty_im_range, false); + TEST_EQUAL (intInt.size(), intMz.size() ) + + TEST_EQUAL (intInt.size(), intIm.size() ) + TEST_EQUAL (intInt.size(), 3) + TEST_REAL_SIMILAR (intInt[0], 0) + TEST_REAL_SIMILAR (intInt[1], 0) + TEST_REAL_SIMILAR (intInt[2], 0) + TEST_REAL_SIMILAR (intMz[0], 101.) // should be middle of window + TEST_REAL_SIMILAR (intMz[1], 103.) // should be middle of window + TEST_REAL_SIMILAR (intMz[2], 105.) // should be middle of window + TEST_REAL_SIMILAR (intIm[0], -1) // should be avg. drift + TEST_REAL_SIMILAR (intIm[1], -1) // should be avg. drift + TEST_REAL_SIMILAR (intIm[2], -1) // should be avg. drift + } + + { + // Test non empty spectrum with no im + std::vector windows, intInt, intMz, intIm; + windows.push_back(101.); + windows.push_back(103.); + windows.push_back(105.); + + DIAHelpers::integrateWindows(spec, windows, 2, intInt, intMz, intIm, empty_im_range); + TEST_EQUAL (intInt.size(), intMz.size() ) + TEST_EQUAL (intInt.size(), intIm.size() ) + TEST_EQUAL (intInt.size(), 3) + TEST_REAL_SIMILAR (intInt[0], 2) + TEST_REAL_SIMILAR (intMz[0], 100.5); + TEST_REAL_SIMILAR (intIm[0], -1); + + // std::cout << "print Int" << std::endl; + // std::copy(intInt.begin(), intInt.end(), + // std::ostream_iterator(std::cout, " ")); + // std::cout << std::endl << "print mz" << intMz.size() << std::endl; + // std::cout << intMz[0] << " " << intMz[1] << " " << intMz[2] << std::endl; + // std::copy(intMz.begin(), intMz.end(), + // std::ostream_iterator(std::cout, " ")); + } + + { + // Test non empty spectrum with ion mobility + std::vector windows, intInt, intMz, intIm; + windows.push_back(102.); + + DIAHelpers::integrateWindows(imSpec, windows, 2, intInt, intMz, intIm, nonempty_im_range); + TEST_EQUAL (intInt.size(), intMz.size() ) + TEST_EQUAL (intInt.size(), intIm.size() ) + TEST_EQUAL (intInt.size(), 1) + TEST_REAL_SIMILAR (intInt[0], 2) + TEST_REAL_SIMILAR (intMz[0], 101.5); + TEST_REAL_SIMILAR (intIm[0], 2.5); + + /* + std::cout << "print Int" << std::endl; + std::copy(intInt.begin(), intInt.end(), + std::ostream_iterator(std::cout, " ")); + std::cout << std::endl << "print mz" << intMz.size() << std::endl; + std::cout << intMz[0] << " " << intMz[1] << " " << intMz[2] << std::endl; + std::copy(intMz.begin(), intMz.end(), + std::ostream_iterator(std::cout, " ")); + */ + } +} +END_SECTION + + +START_SECTION(void integrateWindows(const SpectrumSequence& spectrum, const std::vector & windowsCenter, double width, std::vector & integratedWindowsIntensity, std::vector & integratedWindowsMZ, std::vector & integratedWindowsIm, const RangeMobility& im_range, bool remZero)) +{ + RangeMobility empty_im_range; // empty im range for testing + RangeMobility nonempty_im_range(3.5); + nonempty_im_range.minSpanIfSingular(3); + + { + // Test empty windows + OpenSwath::SpectrumPtr emptySpec(new OpenSwath::Spectrum()); + std::vector specArr; + std::vector windows, intInt, intMz, intIm; + specArr.push_back(emptySpec); + + // message test_exception not working, default to without message for now + //TEST_EXCEPTION_WITH_MESSAGE(Exception::MissingInformation, DIAHelpers::integrateWindows(specArr, windows, 2, intInt, intMz, intIm, empty_im_range), "No windows supplied!"); + TEST_EXCEPTION(Exception::MissingInformation, DIAHelpers::integrateWindows(specArr, windows, 2, intInt, intMz, intIm, empty_im_range)); + } + + + { + // Test empty spectrum (with non empty windows) - remove zeros + OpenSwath::SpectrumPtr emptySpec(new OpenSwath::Spectrum()); + std::vector specArr; + std::vector windows, intInt, intMz, intIm; + specArr.push_back(emptySpec); + + windows.push_back(101.); + windows.push_back(103.); + windows.push_back(105.); + + + DIAHelpers::integrateWindows(specArr, windows, 2, intInt, intMz, intIm, nonempty_im_range, true); + TEST_EQUAL (intInt.empty(), true); + TEST_EQUAL (intIm.empty(), true); + TEST_EQUAL (intMz.empty(), true); + } + + { + // Test empty spectrum (with non empty windows) - Don't remove zeros + OpenSwath::SpectrumPtr emptySpec(new OpenSwath::Spectrum()); + std::vector specArr; + std::vector windows, intInt, intMz, intIm; + specArr.push_back(emptySpec); + + windows.push_back(101.); + windows.push_back(103.); + windows.push_back(105.); + + DIAHelpers::integrateWindows(specArr, windows, 2, intInt, intMz, intIm, empty_im_range, false); + TEST_EQUAL (intInt.size(), intMz.size() ) + TEST_EQUAL (intInt.size(), intIm.size() ) + TEST_EQUAL (intInt.size(), 3) + TEST_REAL_SIMILAR (intInt[0], 0) + TEST_REAL_SIMILAR (intInt[1], 0) + TEST_REAL_SIMILAR (intInt[2], 0) + TEST_REAL_SIMILAR (intMz[0], 101.) // should be middle of window + TEST_REAL_SIMILAR (intMz[1], 103.) // should be middle of window + TEST_REAL_SIMILAR (intMz[2], 105.) // should be middle of window + TEST_REAL_SIMILAR (intIm[0], -1) // should be avg. drift + TEST_REAL_SIMILAR (intIm[1], -1) // should be avg. drift + TEST_REAL_SIMILAR (intIm[2], -1) // should be avg. drift + } + + + { + // Test non empty spectrum with no im + std::vector specArr; + std::vector windows, intInt, intMz, intIm; + windows.push_back(101.); + windows.push_back(103.); + windows.push_back(105.); + + specArr.push_back(spec); + + DIAHelpers::integrateWindows(specArr, windows, 2, intInt, intMz, intIm, empty_im_range); + TEST_EQUAL (intInt.size(), intMz.size() ) + TEST_EQUAL (intInt.size(), intIm.size() ) + TEST_EQUAL (intInt.size(), 3) + TEST_REAL_SIMILAR (intInt[0], 2) + TEST_REAL_SIMILAR (intMz[0], 100.5); + TEST_REAL_SIMILAR (intIm[0], -1); + + // std::cout << "print Int" << std::endl; + // std::copy(intInt.begin(), intInt.end(), + // std::ostream_iterator(std::cout, " ")); + // std::cout << std::endl << "print mz" << intMz.size() << std::endl; + // std::cout << intMz[0] << " " << intMz[1] << " " << intMz[2] << std::endl; + // std::copy(intMz.begin(), intMz.end(), + // std::ostream_iterator(std::cout, " ")); + } + + { + // Test non empty spectrum with ion mobility + std::vector specArr; + std::vector windows, intInt, intMz, intIm; + windows.push_back(102.); + + specArr.push_back(imSpec); + + DIAHelpers::integrateWindows(specArr, windows, 2, intInt, intMz, intIm, nonempty_im_range); + TEST_EQUAL (intInt.size(), intMz.size() ) + TEST_EQUAL (intInt.size(), intIm.size() ) + TEST_EQUAL (intInt.size(), 1) + TEST_REAL_SIMILAR (intInt[0], 2) + TEST_REAL_SIMILAR (intMz[0], 101.5); + TEST_REAL_SIMILAR (intIm[0], 2.5); + + /* + std::cout << "print Int" << std::endl; + std::copy(intInt.begin(), intInt.end(), + std::ostream_iterator(std::cout, " ")); + std::cout << std::endl << "print mz" << intMz.size() << std::endl; + std::cout << intMz[0] << " " << intMz[1] << " " << intMz[2] << std::endl; + std::copy(intMz.begin(), intMz.end(), + std::ostream_iterator(std::cout, " ")); + */ + } } END_SECTION @@ -114,6 +475,15 @@ START_SECTION([EXTRA] void adjustExtractionWindow(double& right, double& left, c TEST_REAL_SIMILAR(left, 500 - 500 * 5 /1e6); TEST_REAL_SIMILAR(right, 500 + 500 * 5 /1e6); } + + // test case where extraction leads to a negative number, should correct this to the 0 bounds (no ppm) + // Note since this is very unlikely with ppm, this functionality is currently not implemented in ppm + { + double left(500.0), right(500.0); + OpenMS::DIAHelpers::adjustExtractionWindow(right, left, 1000, false); + TEST_REAL_SIMILAR(left, 0) + TEST_REAL_SIMILAR(right, 500 + 500) + } } END_SECTION @@ -218,7 +588,6 @@ START_SECTION([EXTRA] simulateSpectrumFromAASequence_test) std::vector masses1; std::vector > tmp, out; OpenMS::DIAHelpers::simulateSpectrumFromAASequence(a, masses1, tmp, &generator); - std::copy(masses1.begin(), masses1.end(), std::ostream_iterator(std::cout, " ")); std::cout << std::endl; diff --git a/src/tests/class_tests/openms/source/DIAPrescoring_test.cpp b/src/tests/class_tests/openms/source/DIAPrescoring_test.cpp index 7d6435f1372..84643ffdeff 100644 --- a/src/tests/class_tests/openms/source/DIAPrescoring_test.cpp +++ b/src/tests/class_tests/openms/source/DIAPrescoring_test.cpp @@ -1,6 +1,6 @@ // Copyright (c) 2002-2023, The OpenMS Team -- EKU Tuebingen, ETH Zurich, and FU Berlin // SPDX-License-Identifier: BSD-3-Clause -// +// // -------------------------------------------------------------------------- // $Maintainer: Timo Sachsenberg $ // $Authors: Witold Wolski $ @@ -11,6 +11,7 @@ #include #include "OpenMS/OPENSWATHALGO/DATAACCESS/MockObjects.h" +#include using namespace std; using namespace OpenMS; @@ -22,6 +23,8 @@ START_TEST(DiaPrescore2, "$Id$") ///////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////// +const std::string ION_MOBILITY_DESCRIPTION = "Ion Mobility"; + DiaPrescore* ptr = nullptr; DiaPrescore* nullPointer = nullptr; @@ -38,7 +41,7 @@ START_SECTION(~DiaPrescore()) } END_SECTION -START_SECTION ( test score function with perfect first transition ) +START_SECTION ( test score function with perfect first transition and ion mobility filtering ) { OpenSwath::LightTransition mock_tr1; mock_tr1.product_mz = 500.; @@ -56,6 +59,7 @@ START_SECTION ( test score function with perfect first transition ) std::vector binaryDataArrayPtrs; OpenSwath::BinaryDataArrayPtr data1(new OpenSwath::BinaryDataArray); OpenSwath::BinaryDataArrayPtr data2(new OpenSwath::BinaryDataArray); + OpenMS::RangeMobility im_range_empty; static const double arr1[] = { 10, 20, 50, 100, 50, 20, 10, // peak at 499 @@ -78,9 +82,10 @@ START_SECTION ( test score function with perfect first transition ) 599.97, 599.98, 599.99, 600.0, 600.01, 600.02, 600.03, 600.97, 600.98, 600.99, 601.0, 601.01, 601.02, 601.03, 601.97, 601.98, 601.99, 602.0, 602.01, 602.02, 602.03, - 602.99, 603.0, 603.01 + 602.99, 603.0, 603.01, }; std::vector mz (arr2, arr2 + sizeof(arr2) / sizeof(double) ); + data1->data = mz; data2->data = intensity; sptr->setMZArray( data1); @@ -92,7 +97,11 @@ START_SECTION ( test score function with perfect first transition ) DiaPrescore diaprescore(0.05); double manhattan = 0., dotprod = 0.; - diaprescore.score(sptr, transitions , dotprod, manhattan); + + std::vector sptrArr; + sptrArr.push_back(sptr); + + diaprescore.score(sptrArr, transitions , im_range_empty, dotprod, manhattan); //std::cout << "dotprod : " << dotprod << std::endl; //std::cout << "manhattan : " << manhattan << std::endl; // >> exp = [240, 74, 39, 15, 0] @@ -121,6 +130,7 @@ START_SECTION ( test score function missing first transition ) mock_tr2.library_intensity = 5.; OpenSwath::SpectrumPtr sptr = (OpenSwath::SpectrumPtr)(new OpenSwath::Spectrum); + std::vector binaryDataArrayPtrs; OpenSwath::BinaryDataArrayPtr data1(new OpenSwath::BinaryDataArray); OpenSwath::BinaryDataArrayPtr data2(new OpenSwath::BinaryDataArray); @@ -162,8 +172,12 @@ START_SECTION ( test score function missing first transition ) transitions.push_back(mock_tr2); DiaPrescore diaprescore(0.05); + OpenMS::RangeMobility im_range_empty; double manhattan = 0., dotprod = 0.; - diaprescore.score(sptr, transitions , dotprod, manhattan); + + std::vector sptrArr; + sptrArr.push_back(sptr); + diaprescore.score(sptrArr, transitions, im_range_empty, dotprod, manhattan); //std::cout << "dotprod : " << dotprod << std::endl; //std::cout << "manhattan : " << manhattan << std::endl; // >> exp = [240, 74, 39, 15, 0] @@ -230,8 +244,13 @@ START_SECTION ( test score function with shifted first transition ) transitions.push_back(mock_tr2); DiaPrescore diaprescore(0.05); + OpenMS::RangeMobility im_range_empty; double manhattan = 0., dotprod = 0.; - diaprescore.score(sptr, transitions , dotprod, manhattan); + + std::vector sptrArr; + sptrArr.push_back(sptr); + + diaprescore.score(sptrArr, transitions, im_range_empty, dotprod, manhattan); //std::cout << "dotprod : " << dotprod << std::endl; //std::cout << "manhattan : " << manhattan << std::endl; // >> exp = [240, 74, 39, 15, 0] @@ -245,7 +264,99 @@ START_SECTION ( test score function with shifted first transition ) } END_SECTION +START_SECTION ( test score function missing first transition due to different ion mobility ) +{ + OpenSwath::LightTransition mock_tr1; + mock_tr1.product_mz = 500.; + mock_tr1.fragment_charge = 1; + mock_tr1.transition_name = "group1"; + mock_tr1.library_intensity = 5.; + + OpenSwath::LightTransition mock_tr2; + mock_tr2.product_mz = 600.; + mock_tr2.fragment_charge = 1; + mock_tr2.transition_name = "group2"; + mock_tr2.library_intensity = 5.; + + double PRECURSOR_ION_MOBILITY = 7; + double ION_MOBILITY_WIDTH = 2; + + OpenSwath::SpectrumPtr sptr = (OpenSwath::SpectrumPtr)(new OpenSwath::Spectrum); + std::vector binaryDataArrayPtrs; + OpenSwath::BinaryDataArrayPtr data1(new OpenSwath::BinaryDataArray); + OpenSwath::BinaryDataArrayPtr data2(new OpenSwath::BinaryDataArray); + OpenSwath::BinaryDataArrayPtr data3(new OpenSwath::BinaryDataArray); + + std::vector intensity { + + 10, 20, 50, 100, 50, 20, 10, // peak at 499 + 3, 7, 15, 30, 15, 7, 3, // peak at 500 + 1, 3, 9, 15, 9, 3, 1, // peak at 501 + 3, 9, 3, // peak at 502 + + 10, 20, 50, 100, 50, 20, 10, // peak at 600 + 3, 7, 15, 30, 15, 7, 3, // peak at 601 + 1, 3, 9, 15, 9, 3, 1, // peak at 602 + 3, 9, 3 // peak at 603 + }; + std::vector mz { + + 498.97, 498.98, 498.99, 499.0, 499.01, 499.02, 499.03, + 499.97, 499.98, 499.99, 500.0, 500.01, 500.02, 500.03, + 500.97, 500.98, 500.99, 501.0, 501.01, 501.02, 501.03, + 501.99, 502.0, 502.01, + + 599.97, 599.98, 599.99, 600.0, 600.01, 600.02, 600.03, + 600.97, 600.98, 600.99, 601.0, 601.01, 601.02, 601.03, + 601.97, 601.98, 601.99, 602.0, 602.01, 602.02, 602.03, + 602.99, 603.0, 603.01 + }; + + std::vector im { + 1, 1, 3, 1, 1, 1, 1, // peak at 499 + 2, 2, 3, 1, 2, 1, 2, // peak at 500 + 1, 2, 1, 2, 3, 2, 2, // peak at 501 + 2, 2, 2, // peak at 502 + + 7, 6, 7, 6, 8, 6, 8, // peak at 600 + 7, 7, 7, 8, 7, 7, 8, // peak at 601 + 8, 6, 8, 8, 6, 6, 6, // peak at 602 + 6, 8, 6 // peak at 603 + }; + + data1->data = mz; + data2->data = intensity; + data3->data = im; + sptr->setMZArray(data1); + sptr->setIntensityArray(data2); + data3->description = ION_MOBILITY_DESCRIPTION; + sptr->getDataArrays().push_back(data3); + + std::vector transitions; + transitions.push_back(mock_tr1); + transitions.push_back(mock_tr2); + + DiaPrescore diaprescore(0.05); + OpenMS::RangeMobility im_range(PRECURSOR_ION_MOBILITY); + im_range.minSpanIfSingular(ION_MOBILITY_WIDTH); + double manhattan = 0., dotprod = 0.; + + std::vector sptrArr; + sptrArr.push_back(sptr); + diaprescore.score(sptrArr, transitions, im_range, dotprod, manhattan); + //std::cout << "dotprod : " << dotprod << std::endl; + //std::cout << "manhattan : " << manhattan << std::endl; + // >> exp = [240, 74, 39, 15, 0] + // >> theo = [1, 0.325757771553019, 0.0678711748364005, 0.0105918703087134, 0.00134955223787482] + // >> from scipy.stats.stats import pearsonr + // >> pearsonr(exp, theo) + // (0.99463189043051314, 0.00047175434098498532) + // + TEST_REAL_SIMILAR(dotprod, 0.627263258948172) + TEST_REAL_SIMILAR(manhattan, 0.984211129641047) +} +END_SECTION + ///////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////// END_TEST - diff --git a/src/tests/class_tests/openms/source/DIAScoring_test.cpp b/src/tests/class_tests/openms/source/DIAScoring_test.cpp index c762acc30f2..8cfd7763e0c 100644 --- a/src/tests/class_tests/openms/source/DIAScoring_test.cpp +++ b/src/tests/class_tests/openms/source/DIAScoring_test.cpp @@ -34,6 +34,67 @@ void getMRMFeatureTest(MockMRMFeature * imrmfeature_test) imrmfeature_test->m_intensity = 1.0; } +OpenSwath::SpectrumPtr prepareIMSpectrum() +{ + OpenSwath::SpectrumPtr sptr = (OpenSwath::SpectrumPtr)(new OpenSwath::Spectrum); + std::vector binaryDataArrayPtrs; + OpenSwath::BinaryDataArrayPtr mz = (OpenSwath::BinaryDataArrayPtr)(new OpenSwath::BinaryDataArray); + OpenSwath::BinaryDataArrayPtr intensity = (OpenSwath::BinaryDataArrayPtr)(new OpenSwath::BinaryDataArray); + OpenSwath::BinaryDataArrayPtr im = (OpenSwath::BinaryDataArrayPtr)(new OpenSwath::BinaryDataArray); + + std::vector intensityVector = { + + 10, 20, 50, 100, 50, 20, 10, // peak at 499 -> 260-20 = 240 intensity within 0.05 Th + 3, 7, 15, 30, 200000, 15, 7, 3, // peak at 500 -> 80-6 = 74 intensity within 0.05 Th + 1, 3, 9, 15, 9, 3, 1, // peak at 501 -> 41-2 = 39 intensity within 0.05 Th + 3, 9, 3, // peak at 502 -> 15 intensity within 0.05 Th + + 10, 20, 50, 100, 50, 20, 10, // peak at 600 -> 260-20 = 240 intensity within 0.05 Th + 3, 7, 15, 30, 15, 7, 3, // peak at 601 -> 80-6 = 74 intensity within 0.05 Th + 1, 3, 9, 15, 9, 3, 1, // peak at 602 -> sum([ 9, 15, 9, 3, 1]) = 37 intensity within 0.05 Th + 3, 9, 200000, 3 // peak at 603 (strong interfering signal (200000) separated out by IM + }; + + std::vector mzVector = { + 498.97, 498.98, 498.99, 499.0, 499.01, 499.02, 499.03, + 499.97, 499.98, 499.99, 500.0, 500.005, 500.01, 500.02, 500.03, // Add interfering 500.025 peak at different IM + 500.97, 500.98, 500.99, 501.0, 501.01, 501.02, 501.03, + 501.99, 502.0, 502.01, + + 599.97, 599.98, 599.99, 600.0, 600.01, 600.02, 600.03, + 600.97, 600.98, 600.99, 601.0, 601.01, 601.02, 601.03, + // note that this peak at 602 is special since it is integrated from + // [(600+2*1.0033548) - 0.025, (600+2*1.0033548) + 0.025] = [601.9817096 to 602.0317096] + 601.97, 601.98, 601.99, 602.0, 602.01, 602.02, 602.03, + 602.99, 603.0, 603.01, 603.01 + }; + + std::vector imVector = { + 2, 3, 2, 3, 4, 2, 3, + 2, 3, 2, 3, 6, 2, 3, 3, // Note the interfering signal has an IM of 6 while all others have IM range of 2-4 + 2, 3, 2, 3, 4, 2, 3, + 2, 3, 3, + + 2, 3, 2, 3, 4, 2, 3, + 2, 3, 2, 3, 4, 2, 3, + 2, 3, 2, 3, 4, 2, 3, + 2, 3, 6, 3, // Note the interfering signal has an IM of 6 while all others have IM range of 2-4 + }; + + TEST_EQUAL(imVector.size(), mzVector.size()); + TEST_EQUAL(imVector.size(), intensityVector.size()); + + mz->data = mzVector; + intensity->data = intensityVector; + im->data = imVector; + im->description = "Ion Mobility"; + + sptr->setMZArray( mz ); + sptr->setIntensityArray( intensity ); + sptr->getDataArrays().push_back( im ); + return sptr; +} + OpenSwath::SpectrumPtr prepareSpectrum() { OpenSwath::SpectrumPtr sptr = (OpenSwath::SpectrumPtr)(new OpenSwath::Spectrum); @@ -41,7 +102,7 @@ OpenSwath::SpectrumPtr prepareSpectrum() OpenSwath::BinaryDataArrayPtr data1 = (OpenSwath::BinaryDataArrayPtr)(new OpenSwath::BinaryDataArray); OpenSwath::BinaryDataArrayPtr data2 = (OpenSwath::BinaryDataArrayPtr)(new OpenSwath::BinaryDataArray); - static const double arr1[] = { + std::vector intensity = { 10, 20, 50, 100, 50, 20, 10, // peak at 499 -> 260-20 = 240 intensity within 0.05 Th 3, 7, 15, 30, 15, 7, 3, // peak at 500 -> 80-6 = 74 intensity within 0.05 Th 1, 3, 9, 15, 9, 3, 1, // peak at 501 -> 41-2 = 39 intensity within 0.05 Th @@ -52,8 +113,7 @@ OpenSwath::SpectrumPtr prepareSpectrum() 1, 3, 9, 15, 9, 3, 1, // peak at 602 -> sum([ 9, 15, 9, 3, 1]) = 37 intensity within 0.05 Th 3, 9, 3 // peak at 603 }; - std::vector intensity (arr1, arr1 + sizeof(arr1) / sizeof(arr1[0]) ); - static const double arr2[] = { + std::vector mz = { 498.97, 498.98, 498.99, 499.0, 499.01, 499.02, 499.03, 499.97, 499.98, 499.99, 500.0, 500.01, 500.02, 500.03, 500.97, 500.98, 500.99, 501.0, 501.01, 501.02, 501.03, @@ -61,21 +121,36 @@ OpenSwath::SpectrumPtr prepareSpectrum() 599.97, 599.98, 599.99, 600.0, 600.01, 600.02, 600.03, 600.97, 600.98, 600.99, 601.0, 601.01, 601.02, 601.03, - // note that this peak at 602 is special since it is integrated from + // note that this peak at 602 is special since it is integrated from // [(600+2*1.0033548) - 0.025, (600+2*1.0033548) + 0.025] = [601.9817096 to 602.0317096] 601.97, 601.98, 601.99, 602.0, 602.01, 602.02, 602.03, 602.99, 603.0, 603.01 }; - std::vector mz (arr2, arr2 + sizeof(arr2) / sizeof(arr2[0]) ); data1->data = mz; data2->data = intensity; - sptr->setMZArray(data1); + sptr->setMZArray( data1 ); sptr->setIntensityArray( data2 ); return sptr; } -OpenSwath::SpectrumPtr prepareShiftedSpectrum() +SpectrumSequence prepareSpectrumSequence() +{ + OpenSwath::SpectrumPtr sptr = prepareSpectrum(); + SpectrumSequence out; + out.push_back(sptr); + return out; +} + +SpectrumSequence prepareIMSpectrumSequence() +{ + OpenSwath::SpectrumPtr sptr = prepareIMSpectrum(); + SpectrumSequence out; + out.push_back(sptr); + return out; +} + +SpectrumSequence prepareShiftedSpectrum() { OpenSwath::SpectrumPtr sptr = prepareSpectrum(); // shift the peaks by a fixed amount in ppm @@ -87,9 +162,29 @@ OpenSwath::SpectrumPtr prepareShiftedSpectrum() { sptr->getMZArray()->data[i] += sptr->getMZArray()->data[i] / 1000000 * 10; // shift second peak by 10 ppm } - return sptr; + SpectrumSequence out; + out.push_back(sptr); + return out; +} + +SpectrumSequence prepareShiftedSpectrumIM() +{ + OpenSwath::SpectrumPtr sptr = prepareIMSpectrum(); + // shift the peaks by a fixed amount in ppm + for (std::size_t i = 0; i < sptr->getMZArray()->data.size() / 2.0; i++) + { + sptr->getMZArray()->data[i] += sptr->getMZArray()->data[i] / 1000000 * 15; // shift first peak by 15 ppm + } + for (std::size_t i = sptr->getMZArray()->data.size() / 2.0; i < sptr->getMZArray()->data.size(); i++) + { + sptr->getMZArray()->data[i] += sptr->getMZArray()->data[i] / 1000000 * 10; // shift second peak by 10 ppm + } + SpectrumSequence out; + out.push_back(sptr); + return out; } + START_TEST(DIAScoring, "$Id$") ///////////////////////////////////////////////////////////// @@ -196,21 +291,23 @@ mock_tr2.product_mz = 600; mock_tr2.fragment_charge = 1; mock_tr2.transition_name = "group2"; -START_SECTION([EXTRA] forward void dia_isotope_scores(const std::vector & transitions, SpectrumType spectrum, OpenSwath::IMRMFeature * mrmfeature, int putative_fragment_charge, double & isotope_corr, double & isotope_overlap)) +START_SECTION([EXTRA] forward void dia_isotope_scores(const std::vector & transitions, std::vector spectrum, OpenSwath::IMRMFeature * mrmfeature, int putative_fragment_charge, const RangeMobility& im_range, double & isotope_corr, double & isotope_overlap)) { - OpenSwath::SpectrumPtr sptr = prepareSpectrum(); + SpectrumSequence sptr = prepareSpectrumSequence(); MockMRMFeature * imrmfeature_test = new MockMRMFeature(); getMRMFeatureTest(imrmfeature_test); imrmfeature_test->m_intensity = 0.7f; std::vector transitions; - // Try with transition at 600 m/z + OpenMS::RangeMobility empty_im_range; + + // Try with transition at 600 m/z transitions.push_back(mock_tr2); DIAScoring diascoring; diascoring.setParameters(p_dia); double isotope_corr = 0, isotope_overlap = 0; - diascoring.dia_isotope_scores(transitions, sptr, imrmfeature_test, isotope_corr, isotope_overlap); + diascoring.dia_isotope_scores(transitions, sptr, imrmfeature_test, empty_im_range, isotope_corr, isotope_overlap); // >> exp = [240, 74, 37, 15, 0] // >> theo = [1, 0.325757771553019, 0.0678711748364005, 0.0105918703087134, 0.00134955223787482] @@ -224,9 +321,40 @@ START_SECTION([EXTRA] forward void dia_isotope_scores(const std::vector & transitions, SpectrumType spectrum, OpenSwath::IMRMFeature * mrmfeature, int putative_fragment_charge, double & isotope_corr, double & isotope_overlap)) + +START_SECTION([EXTRA] forward with IM void dia_isotope_scores(const std::vector & transitions, std::vector spectrum, OpenSwath::IMRMFeature * mrmfeature, int putative_fragment_charge, const RangeMobility& im_range, double & isotope_corr, double & isotope_overlap)) { - OpenSwath::SpectrumPtr sptr = prepareSpectrum(); + SpectrumSequence sptr = prepareIMSpectrumSequence(); + MockMRMFeature * imrmfeature_test = new MockMRMFeature(); + getMRMFeatureTest(imrmfeature_test); + imrmfeature_test->m_intensity = 0.7f; + std::vector transitions; + + // im_range 2-4 + OpenMS::RangeMobility im_range(3); + im_range.minSpanIfSingular(2); // im_range from 2-4 + + // Try with transition at 600 m/z + transitions.push_back(mock_tr2); + + DIAScoring diascoring; + diascoring.setParameters(p_dia); + double isotope_corr = 0, isotope_overlap = 0; + diascoring.dia_isotope_scores(transitions, sptr, imrmfeature_test, im_range, isotope_corr, isotope_overlap); + + // Should be same as above because the IM is filtered out + TEST_REAL_SIMILAR(isotope_corr, 0.995335798317618) + TEST_REAL_SIMILAR(isotope_overlap, 0.0) + delete imrmfeature_test; +} +END_SECTION + + + +START_SECTION([EXTRA] forward negative charge: void dia_isotope_scores(const std::vector & transitions, std::vector spectrum, OpenSwath::IMRMFeature * mrmfeature, int putative_fragment_charge, const RangeMobility& im_range, double & isotope_corr, double & isotope_overlap)) +{ + RangeMobility empty_im_range; + SpectrumSequence sptr = prepareSpectrumSequence(); MockMRMFeature * imrmfeature_test = new MockMRMFeature(); getMRMFeatureTest(imrmfeature_test); imrmfeature_test->m_intensity = 0.7f; @@ -239,7 +367,7 @@ START_SECTION([EXTRA] forward negative charge: void dia_isotope_scores(const std DIAScoring diascoring; diascoring.setParameters(p_dia); double isotope_corr = 0, isotope_overlap = 0; - diascoring.dia_isotope_scores(transitions, sptr, imrmfeature_test, isotope_corr, isotope_overlap); + diascoring.dia_isotope_scores(transitions, sptr, imrmfeature_test, empty_im_range, isotope_corr, isotope_overlap); // >> exp = [240, 74, 37, 15, 0] // >> theo = [1, 0.325757771553019, 0.0678711748364005, 0.0105918703087134, 0.00134955223787482] @@ -253,9 +381,10 @@ START_SECTION([EXTRA] forward negative charge: void dia_isotope_scores(const std } END_SECTION -START_SECTION([EXTRA] forward negative charge: void dia_isotope_scores(const std::vector & transitions, SpectrumType spectrum, OpenSwath::IMRMFeature * mrmfeature, int putative_fragment_charge, double & isotope_corr, double & isotope_overlap)) +START_SECTION([EXTRA] forward negative charge: void dia_isotope_scores(const std::vector & transitions, std::vector spectrum, OpenSwath::IMRMFeature * mrmfeature, int putative_fragment_charge, const RangeMobility& im_range, double & isotope_corr, double & isotope_overlap)) { - OpenSwath::SpectrumPtr sptr = prepareSpectrum(); + OpenMS::RangeMobility empty_im_range; + SpectrumSequence sptr = prepareSpectrumSequence(); MockMRMFeature * imrmfeature_test = new MockMRMFeature(); getMRMFeatureTest(imrmfeature_test); imrmfeature_test->m_intensity = 0.7f; @@ -268,7 +397,7 @@ START_SECTION([EXTRA] forward negative charge: void dia_isotope_scores(const std DIAScoring diascoring; diascoring.setParameters(p_dia); double isotope_corr = 0, isotope_overlap = 0; - diascoring.dia_isotope_scores(transitions, sptr, imrmfeature_test, isotope_corr, isotope_overlap); + diascoring.dia_isotope_scores(transitions, sptr, imrmfeature_test, empty_im_range, isotope_corr, isotope_overlap); // >> exp = [240, 74, 37, 15, 0] // >> theo = [1, 0.325757771553019, 0.0678711748364005, 0.0105918703087134, 0.00134955223787482] @@ -282,22 +411,23 @@ START_SECTION([EXTRA] forward negative charge: void dia_isotope_scores(const std } END_SECTION -START_SECTION([EXTRA] backward void dia_isotope_scores(const std::vector & transitions, SpectrumType spectrum, OpenSwath::IMRMFeature * mrmfeature, int putative_fragment_charge, double & isotope_corr, double & isotope_overlap)) +START_SECTION([EXTRA] backward void dia_isotope_scores(const std::vector & transitions, std::vector spectrum, OpenSwath::IMRMFeature * mrmfeature, int putative_fragment_charge, const RangeMobility& im_range, double & isotope_corr, double & isotope_overlap)) { - OpenSwath::SpectrumPtr sptr = prepareSpectrum(); + OpenMS::RangeMobility empty_im_range; + SpectrumSequence sptr = prepareSpectrumSequence(); MockMRMFeature * imrmfeature_test = new MockMRMFeature(); getMRMFeatureTest(imrmfeature_test); imrmfeature_test->m_intensity = 0.3f; std::vector transitions; - // Try with transition at 500 m/z + // Try with transition at 500 m/z // This peak is not monoisotopic (e.g. at 499 there is another, more intense, peak) transitions.push_back(mock_tr1); DIAScoring diascoring; diascoring.setParameters(p_dia); double isotope_corr = 0, isotope_overlap = 0; - diascoring.dia_isotope_scores(transitions, sptr, imrmfeature_test, isotope_corr, isotope_overlap); + diascoring.dia_isotope_scores(transitions, sptr, imrmfeature_test, empty_im_range, isotope_corr, isotope_overlap); // >> exp = [74, 39, 15, 0, 0] // >> theo = [1, 0.266799519434277, 0.0486475002325161, 0.0066525896497495, 0.000747236543377621] @@ -310,9 +440,69 @@ START_SECTION([EXTRA] backward void dia_isotope_scores(const std::vector &transitions, SpectrumType spectrum, OpenSwath::IMRMFeature *mrmfeature, double &isotope_corr, double &isotope_overlap) ) +START_SECTION([EXTRA] backward with IM void dia_isotope_scores(const std::vector & transitions, std::vector spectrum, OpenSwath::IMRMFeature * mrmfeature, int putative_fragment_charge, const RangeMobility& im_range, double & isotope_corr, double & isotope_overlap)) { - OpenSwath::SpectrumPtr sptr = prepareSpectrum(); + // IM range from 2-4 + OpenMS::RangeMobility im_range(3); + im_range.minSpanIfSingular(2); + + SpectrumSequence sptr = prepareIMSpectrumSequence(); + MockMRMFeature * imrmfeature_test = new MockMRMFeature(); + getMRMFeatureTest(imrmfeature_test); + imrmfeature_test->m_intensity = 0.3f; + std::vector transitions; + + // Try with transition at 500 m/z + // This peak is not monoisotopic (e.g. at 499 there is another, more intense, peak) + transitions.push_back(mock_tr1); + + DIAScoring diascoring; + diascoring.setParameters(p_dia); + double isotope_corr = 0, isotope_overlap = 0; + diascoring.dia_isotope_scores(transitions, sptr, imrmfeature_test, im_range, isotope_corr, isotope_overlap); + + // >> exp = [74, 39, 15, 0, 0] + // >> theo = [1, 0.266799519434277, 0.0486475002325161, 0.0066525896497495, 0.000747236543377621] + // >> from scipy.stats.stats import pearsonr + // >> pearsonr(exp, theo) + // (0.959570883150479, 0.0096989307464742554) + TEST_REAL_SIMILAR(isotope_corr, 0.959692139694113) + TEST_REAL_SIMILAR(isotope_overlap, 1.0) + delete imrmfeature_test; +} +END_SECTION + + +START_SECTION ( void dia_isotope_scores(const std::vector< TransitionType > &transitions, std::vector spectrum, OpenSwath::IMRMFeature *mrmfeature, const RangeMobility& im_range, double &isotope_corr, double &isotope_overlap) ) +{ + OpenMS::RangeMobility empty_im_range; + SpectrumSequence sptr = prepareSpectrumSequence(); + + MockMRMFeature * imrmfeature_test = new MockMRMFeature(); + getMRMFeatureTest(imrmfeature_test); + + // create transitions, e.g. library intensity + std::vector transitions; + transitions.push_back(mock_tr1); + transitions.push_back(mock_tr2); + + DIAScoring diascoring; + diascoring.setParameters(p_dia); + double isotope_corr = 0, isotope_overlap = 0; + diascoring.dia_isotope_scores(transitions, sptr, imrmfeature_test, empty_im_range, isotope_corr, isotope_overlap); + + // see above for the two individual numbers (forward and backward) + TEST_REAL_SIMILAR(isotope_corr, 0.995335798317618 * 0.7 + 0.959692139694113 * 0.3) + TEST_REAL_SIMILAR(isotope_overlap, 0.0 * 0.7 + 1.0 * 0.3) + delete imrmfeature_test; +} +END_SECTION + +START_SECTION ( [EXTRA] with IM void dia_isotope_scores(const std::vector< TransitionType > &transitions, std::vector spectrum, OpenSwath::IMRMFeature *mrmfeature, const RangeMobility& im_range, double &isotope_corr, double &isotope_overlap) ) +{ + OpenMS::RangeMobility im_range(3); + im_range.minSpanIfSingular(2); + SpectrumSequence sptr = prepareIMSpectrumSequence(); MockMRMFeature * imrmfeature_test = new MockMRMFeature(); getMRMFeatureTest(imrmfeature_test); @@ -325,7 +515,7 @@ START_SECTION ( void dia_isotope_scores(const std::vector< TransitionType > &tra DIAScoring diascoring; diascoring.setParameters(p_dia); double isotope_corr = 0, isotope_overlap = 0; - diascoring.dia_isotope_scores(transitions, sptr, imrmfeature_test, isotope_corr, isotope_overlap); + diascoring.dia_isotope_scores(transitions, sptr, imrmfeature_test, im_range, isotope_corr, isotope_overlap); // see above for the two individual numbers (forward and backward) TEST_REAL_SIMILAR(isotope_corr, 0.995335798317618 * 0.7 + 0.959692139694113 * 0.3) @@ -334,14 +524,16 @@ START_SECTION ( void dia_isotope_scores(const std::vector< TransitionType > &tra } END_SECTION -START_SECTION(void dia_ms1_isotope_scores(double precursor_mz, SpectrumPtrType spectrum, size_t charge_state, +START_SECTION(void dia_ms1_isotope_scores(double precursor_mz, std::vector spectrum, size_t charge_state, const RangeMobility& im_range, double& isotope_corr, double& isotope_overlap)) { - OpenSwath::SpectrumPtr sptr = prepareSpectrum(); + SpectrumSequence sptr = prepareSpectrumSequence(); DIAScoring diascoring; diascoring.setParameters(p_dia); + OpenMS::RangeMobility empty_im_range; + // Check for charge 1+ and m/z at 500 { size_t precursor_charge_state = 1; @@ -349,7 +541,7 @@ START_SECTION(void dia_ms1_isotope_scores(double precursor_mz, SpectrumPtrType s double isotope_corr = 0, isotope_overlap = 0; diascoring - .dia_ms1_isotope_scores_averagine(precursor_mz, sptr, isotope_corr, isotope_overlap, precursor_charge_state); + .dia_ms1_isotope_scores_averagine(precursor_mz, sptr, precursor_charge_state, empty_im_range, isotope_corr, isotope_overlap); // see above for the two individual numbers (forward and backward) TEST_REAL_SIMILAR(isotope_corr, 0.959692139694113) @@ -363,7 +555,7 @@ START_SECTION(void dia_ms1_isotope_scores(double precursor_mz, SpectrumPtrType s double isotope_corr = 0, isotope_overlap = 0; diascoring - .dia_ms1_isotope_scores_averagine(precursor_mz, sptr, isotope_corr, isotope_overlap, precursor_charge_state); + .dia_ms1_isotope_scores_averagine(precursor_mz, sptr, precursor_charge_state, empty_im_range, isotope_corr, isotope_overlap); // >>> theo = [0.57277789564886, 0.305415548811564, 0.0952064968352544, 0.0218253361702587, 0.00404081869309618] // >>> exp = [74, 0, 39, 0, 15] @@ -380,7 +572,7 @@ START_SECTION(void dia_ms1_isotope_scores(double precursor_mz, SpectrumPtrType s double isotope_corr = 0, isotope_overlap = 0; diascoring - .dia_ms1_isotope_scores_averagine(precursor_mz, sptr, isotope_corr, isotope_overlap, precursor_charge_state); + .dia_ms1_isotope_scores_averagine(precursor_mz, sptr, precursor_charge_state, empty_im_range, isotope_corr, isotope_overlap); // >> exp = [240, 74, 39, 15, 0] // >> theo = [0.755900817146293, 0.201673974754608, 0.0367726851778834, 0.00502869795238462, 0.000564836713740715] @@ -393,9 +585,108 @@ START_SECTION(void dia_ms1_isotope_scores(double precursor_mz, SpectrumPtrType s } END_SECTION -START_SECTION (void dia_massdiff_score(const std::vector< TransitionType > &transitions, SpectrumType spectrum, const std::vector< double > &normalized_library_intensity, double &ppm_score, double &ppm_score_weighted) ) +START_SECTION([EXTRA] with IM void dia_ms1_isotope_scores(double precursor_mz, std::vector spectrum, size_t charge_state, const RangeMobility& im_range, + double& isotope_corr, double& isotope_overlap)) { - OpenSwath::SpectrumPtr sptr = prepareShiftedSpectrum(); + SpectrumSequence sptr = prepareIMSpectrumSequence(); + + DIAScoring diascoring; + diascoring.setParameters(p_dia); + + // IM range 2-4 + OpenMS::RangeMobility im_range(3); + im_range.minSpanIfSingular(2); + + // Check for charge 1+ and m/z at 500 + { + size_t precursor_charge_state = 1; + double precursor_mz = 500; + + double isotope_corr = 0, isotope_overlap = 0; + diascoring + .dia_ms1_isotope_scores_averagine(precursor_mz, sptr, precursor_charge_state, im_range, isotope_corr, isotope_overlap); + + // see above for the two individual numbers (forward and backward) + TEST_REAL_SIMILAR(isotope_corr, 0.959692139694113) + TEST_REAL_SIMILAR(isotope_overlap, 240/74.0) + } + + // Check if charge state is assumed 2+ + { + size_t precursor_charge_state = 2; + double precursor_mz = 500; + + double isotope_corr = 0, isotope_overlap = 0; + diascoring + .dia_ms1_isotope_scores_averagine(precursor_mz, sptr, precursor_charge_state, im_range, isotope_corr, isotope_overlap); + + // >>> theo = [0.57277789564886, 0.305415548811564, 0.0952064968352544, 0.0218253361702587, 0.00404081869309618] + // >>> exp = [74, 0, 39, 0, 15] + // >>> pearsonr(exp, theo) + // (0.68135233883093205, 0.20528953804781694) + TEST_REAL_SIMILAR(isotope_corr, 0.680028918385546) + TEST_REAL_SIMILAR(isotope_overlap, 240/74.0) + } + + // Check and confirm that monoisotopic is at m/z 499 + { + size_t precursor_charge_state = 1; + double precursor_mz = 499; + + double isotope_corr = 0, isotope_overlap = 0; + diascoring + .dia_ms1_isotope_scores_averagine(precursor_mz, sptr, precursor_charge_state, im_range, isotope_corr, isotope_overlap); + + // >> exp = [240, 74, 39, 15, 0] + // >> theo = [0.755900817146293, 0.201673974754608, 0.0367726851778834, 0.00502869795238462, 0.000564836713740715] + // >> from scipy.stats.stats import pearsonr + // >> pearsonr(exp, theo) + // (0.99463189043051314, 0.00047175434098498532) + TEST_REAL_SIMILAR(isotope_corr, 0.995485552148335) + TEST_REAL_SIMILAR(isotope_overlap, 0.0) // monoisotopic + } +} +END_SECTION + + + +START_SECTION (void dia_massdiff_score(const std::vector< TransitionType > &transitions, const SpectrumSequence& spectrum, const std::vector< double > &normalized_library_intensity, const RangeMobility& im_range, double &ppm_score, double &ppm_score_weighted) ) +{ + SpectrumSequence sptr = prepareShiftedSpectrum(); + + OpenMS::RangeMobility empty_im_range; + MockMRMFeature * imrmfeature_test = new MockMRMFeature(); + getMRMFeatureTest(imrmfeature_test); + delete imrmfeature_test; + + // create transitions, e.g. library intensity + std::vector transitions; + transitions.push_back(mock_tr1); + transitions.push_back(mock_tr2); + + DIAScoring diascoring; + diascoring.setParameters(p_dia_large); + double ppm_score = 0, ppm_score_weighted = 0; + std::vector normalized_library_intensity; + normalized_library_intensity.push_back(0.7); + normalized_library_intensity.push_back(0.3); + std::vector ppm_errors; + diascoring.dia_massdiff_score(transitions, sptr, normalized_library_intensity, empty_im_range, ppm_score, ppm_score_weighted, ppm_errors); + + TEST_REAL_SIMILAR(ppm_score, (15 + 10) / 2.0); // 15 ppm and 10 ppm + TEST_REAL_SIMILAR(ppm_score_weighted, 15 * 0.7 + 10* 0.3); // weighted +} +END_SECTION + + +START_SECTION ([EXTRA] with IM void dia_massdiff_score(const std::vector< TransitionType > &transitions, const SpectrumSequence& spectrum, const std::vector< double > &normalized_library_intensity, const RangeMobility& im_range, double &ppm_score, double &ppm_score_weighted) ) +{ + SpectrumSequence sptr = prepareShiftedSpectrumIM(); + + // IM range from 2-4 + OpenMS::RangeMobility im_range(3); + im_range.minSpanIfSingular(2); + MockMRMFeature * imrmfeature_test = new MockMRMFeature(); getMRMFeatureTest(imrmfeature_test); @@ -413,38 +704,64 @@ START_SECTION (void dia_massdiff_score(const std::vector< TransitionType > &tran normalized_library_intensity.push_back(0.7); normalized_library_intensity.push_back(0.3); std::vector ppm_errors; - diascoring.dia_massdiff_score(transitions, sptr, normalized_library_intensity, ppm_score, ppm_score_weighted, ppm_errors); + diascoring.dia_massdiff_score(transitions, sptr, normalized_library_intensity, im_range, ppm_score, ppm_score_weighted, ppm_errors); TEST_REAL_SIMILAR(ppm_score, (15 + 10) / 2.0); // 15 ppm and 10 ppm TEST_REAL_SIMILAR(ppm_score_weighted, 15 * 0.7 + 10* 0.3); // weighted } END_SECTION -START_SECTION ( bool DIAScoring::dia_ms1_massdiff_score(double precursor_mz, transitions, SpectrumType spectrum, double& ppm_score) ) -{ - OpenSwath::SpectrumPtr sptr = prepareShiftedSpectrum(); + +START_SECTION ( bool DIAScoring::dia_ms1_massdiff_score(double precursor_mz, transitions, const SpectrumSequence& spectrum, RangeMobility& im_range, double& ppm_score) ) +{ + SpectrumSequence sptr = prepareShiftedSpectrum(); + DIAScoring diascoring; + OpenMS::RangeMobility empty_imRange; + diascoring.setParameters(p_dia_large); + double ppm_score = 0; + + TEST_EQUAL(diascoring.dia_ms1_massdiff_score(500.0, sptr, empty_imRange, ppm_score), true); + TEST_REAL_SIMILAR(ppm_score, 15); // 15 ppm shifted + + TEST_EQUAL(diascoring.dia_ms1_massdiff_score(600.0, sptr, empty_imRange, ppm_score), true); + TEST_REAL_SIMILAR(ppm_score, 10); // 10 ppm shifted + + TEST_EQUAL(diascoring.dia_ms1_massdiff_score(100.0, sptr, empty_imRange, ppm_score), false); + TEST_REAL_SIMILAR(ppm_score, 0.5 * 1000000 / 100.0); // not present +} +END_SECTION + +START_SECTION ( [EXTRA] with IM bool DIAScoring::dia_ms1_massdiff_score(double precursor_mz, transitions, const SpectrumSequence& spectrum, RangeMobility& im_range, double& ppm_score) ) +{ + SpectrumSequence sptr = prepareShiftedSpectrumIM(); DIAScoring diascoring; + + // im_range 2-4 + OpenMS::RangeMobility imRange(3); + imRange.minSpanIfSingular(2); diascoring.setParameters(p_dia_large); double ppm_score = 0; - TEST_EQUAL(diascoring.dia_ms1_massdiff_score(500.0, sptr, ppm_score), true); + TEST_EQUAL(diascoring.dia_ms1_massdiff_score(500.0, sptr, imRange, ppm_score), true); TEST_REAL_SIMILAR(ppm_score, 15); // 15 ppm shifted - TEST_EQUAL(diascoring.dia_ms1_massdiff_score(600.0, sptr, ppm_score), true); + TEST_EQUAL(diascoring.dia_ms1_massdiff_score(600.0, sptr, imRange, ppm_score), true); TEST_REAL_SIMILAR(ppm_score, 10); // 10 ppm shifted - TEST_EQUAL(diascoring.dia_ms1_massdiff_score(100.0, sptr, ppm_score), false); + TEST_EQUAL(diascoring.dia_ms1_massdiff_score(100.0, sptr, imRange, ppm_score), false); TEST_REAL_SIMILAR(ppm_score, 0.5 * 1000000 / 100.0); // not present } END_SECTION -START_SECTION ( void dia_by_ion_score(SpectrumType spectrum, AASequence &sequence, int charge, double &bseries_score, double &yseries_score) ) + +START_SECTION ( void dia_by_ion_score(std::vector spectrum, AASequence &sequence, int charge, RangeMobility& im_range, double &bseries_score, double &yseries_score) ) { OpenSwath::SpectrumPtr sptr = (OpenSwath::SpectrumPtr)(new OpenSwath::Spectrum); std::vector binaryDataArrayPtrs; OpenSwath::BinaryDataArrayPtr data1 = (OpenSwath::BinaryDataArrayPtr)(new OpenSwath::BinaryDataArray); OpenSwath::BinaryDataArrayPtr data2 = (OpenSwath::BinaryDataArrayPtr)(new OpenSwath::BinaryDataArray); + OpenMS::RangeMobility empty_imRange; std::vector intensity(6, 100); std::vector mz { // four of the naked b/y ions @@ -470,7 +787,9 @@ START_SECTION ( void dia_by_ion_score(SpectrumType spectrum, AASequence &sequenc AASequence a = AASequence::fromString(sequence); double bseries_score = 0, yseries_score = 0; - diascoring.dia_by_ion_score(sptr, a, 1, bseries_score, yseries_score); + SpectrumSequence sptrArr; + sptrArr.push_back(sptr); + diascoring.dia_by_ion_score(sptrArr, a, 1, empty_imRange, bseries_score, yseries_score); TEST_REAL_SIMILAR (bseries_score, 2); TEST_REAL_SIMILAR (yseries_score, 2); @@ -478,14 +797,14 @@ START_SECTION ( void dia_by_ion_score(SpectrumType spectrum, AASequence &sequenc // now add a modification to the sequence a.setModification(1, "Phospho" ); // modify the Y bseries_score = 0, yseries_score = 0; - diascoring.dia_by_ion_score(sptr, a, 1, bseries_score, yseries_score); + diascoring.dia_by_ion_score(sptrArr, a, 1, empty_imRange, bseries_score, yseries_score); TEST_REAL_SIMILAR (bseries_score, 1); TEST_REAL_SIMILAR (yseries_score, 3); } END_SECTION -START_SECTION( void score_with_isotopes(SpectrumType spectrum, const std::vector< TransitionType > &transitions, double &dotprod, double &manhattan)) +START_SECTION( void score_with_isotopes(std::vector spectrum, const std::vector< TransitionType > &transitions, double &dotprod, double &manhattan)) { OpenSwath::LightTransition mock_tr1; mock_tr1.product_mz = 500.; @@ -499,7 +818,39 @@ START_SECTION( void score_with_isotopes(SpectrumType spectrum, const std::vector mock_tr2.transition_name = "group2"; mock_tr2.library_intensity = 5.; - OpenSwath::SpectrumPtr sptr = prepareSpectrum(); + SpectrumSequence sptr = prepareSpectrumSequence(); + + std::vector transitions; + transitions.push_back(mock_tr1); + transitions.push_back(mock_tr2); + + DIAScoring diascoring; + diascoring.setParameters(p_dia); + OpenMS::RangeMobility empty_imRange; + + double dotprod, manhattan; + diascoring.score_with_isotopes(sptr,transitions, empty_imRange, dotprod,manhattan); + + TEST_REAL_SIMILAR (dotprod, 0.43738312458795); + TEST_REAL_SIMILAR (manhattan, 0.55743322213368); +} +END_SECTION + +START_SECTION( [EXTRA] with IM void score_with_isotopes(std::vector spectrum, const std::vector< TransitionType > &transitions, double &dotprod, double &manhattan)) +{ + OpenSwath::LightTransition mock_tr1; + mock_tr1.product_mz = 500.; + mock_tr1.fragment_charge = 1; + mock_tr1.transition_name = "group1"; + mock_tr1.library_intensity = 5.; + + OpenSwath::LightTransition mock_tr2; + mock_tr2.product_mz = 600.; + mock_tr2.fragment_charge = 1; + mock_tr2.transition_name = "group2"; + mock_tr2.library_intensity = 5.; + + SpectrumSequence sptr = prepareIMSpectrumSequence(); std::vector transitions; transitions.push_back(mock_tr1); @@ -508,8 +859,12 @@ START_SECTION( void score_with_isotopes(SpectrumType spectrum, const std::vector DIAScoring diascoring; diascoring.setParameters(p_dia); + // im range from 2-4 + OpenMS::RangeMobility imRange(3); + imRange.minSpanIfSingular(2); + double dotprod, manhattan; - diascoring.score_with_isotopes(sptr,transitions,dotprod,manhattan); + diascoring.score_with_isotopes(sptr,transitions, imRange, dotprod,manhattan); TEST_REAL_SIMILAR (dotprod, 0.43738312458795); TEST_REAL_SIMILAR (manhattan, 0.55743322213368); diff --git a/src/tests/class_tests/openms/source/IonMobilityScoring_test.cpp b/src/tests/class_tests/openms/source/IonMobilityScoring_test.cpp index c9481bf6bac..85114f36053 100644 --- a/src/tests/class_tests/openms/source/IonMobilityScoring_test.cpp +++ b/src/tests/class_tests/openms/source/IonMobilityScoring_test.cpp @@ -190,12 +190,11 @@ OpenSwath::SpectrumPtr ms1spec(new OpenSwath::Spectrum()); } START_SECTION(([EXTRA] - static void driftScoring(OpenSwath::SpectrumPtr spectrum, + static void driftScoring(const SpectrumSequence& spectrum, const std::vector & transitions, OpenSwath_Scores & scores, - const double drift_lower, - const double drift_upper, const double drift_target, + const RangeMobility im_range, const double dia_extraction_window_, const bool dia_extraction_ppm_, const bool use_spline, @@ -203,40 +202,51 @@ START_SECTION(([EXTRA] { OpenSwath_Scores scores; - double drift_lower = 0.5; - double drift_upper = 1.5; double drift_target = 1.0; + OpenMS::RangeMobility im_range_1(1.0); + im_range_1.minSpanIfSingular(1.); double im_drift_extra_pcnt_ = 0.25; double dia_extract_window_ = 0.3; bool dia_extraction_ppm_ = false; - OpenSwath::SpectrumPtr drift_spectrum(new OpenSwath::Spectrum()); - OpenSwath::BinaryDataArrayPtr ion_mobility(new OpenSwath::BinaryDataArray); - drift_spectrum->getDataArrays().push_back( ion_mobility ); - - IonMobilityScoring::driftScoring(drift_spectrum, transitions, scores, - drift_lower, drift_upper, drift_target, - dia_extract_window_, dia_extraction_ppm_, - false, im_drift_extra_pcnt_); + // Test #1: Empty Spectrum + OpenSwath::SpectrumPtr drift_spectrum(new OpenSwath::Spectrum()); OpenSwath::BinaryDataArrayPtr mass(new OpenSwath::BinaryDataArray); OpenSwath::BinaryDataArrayPtr intensity(new OpenSwath::BinaryDataArray); + OpenSwath::BinaryDataArrayPtr ion_mobility(new OpenSwath::BinaryDataArray); + drift_spectrum->setMZArray( mass); drift_spectrum->setIntensityArray( intensity); + drift_spectrum->getDataArrays().push_back( ion_mobility ); - IonMobilityScoring::driftScoring(drift_spectrum, transitions, scores, - drift_lower, drift_upper, drift_target, + std::vector sptrArr; + sptrArr.push_back(drift_spectrum); + + IonMobilityScoring::driftScoring(sptrArr, transitions, scores, + drift_target, im_range_1, dia_extract_window_, dia_extraction_ppm_, false, im_drift_extra_pcnt_); + TEST_REAL_SIMILAR(scores.im_drift, 0); + TEST_REAL_SIMILAR(scores.im_drift_weighted, 0); + TEST_REAL_SIMILAR(scores.im_delta_score, 0); + TEST_REAL_SIMILAR(scores.im_xcorr_shape_score, 0) + TEST_REAL_SIMILAR(scores.im_xcorr_coelution_score, 0) + + // Test #2: IM Scores (Condition 1/2) drift_spectrum = spec; + std::vector sptrArr2; + sptrArr2.push_back(spec); + // Test integrity of spectrum TEST_EQUAL(drift_spectrum->getMZArray()->data.size(), 24) TEST_EQUAL(drift_spectrum->getMZArray()->data.size(), drift_spectrum->getIntensityArray()->data.size()) TEST_EQUAL(drift_spectrum->getMZArray()->data.size(), drift_spectrum->getDriftTimeArray()->data.size()) - IonMobilityScoring::driftScoring(drift_spectrum, transitions, scores, - drift_lower, drift_upper, drift_target, + + IonMobilityScoring::driftScoring(sptrArr2, transitions, scores, + drift_target, im_range_1, dia_extract_window_, dia_extraction_ppm_, false, im_drift_extra_pcnt_); @@ -247,9 +257,10 @@ START_SECTION(([EXTRA] TEST_REAL_SIMILAR(scores.im_xcorr_shape_score, 0.892124778448826) TEST_REAL_SIMILAR(scores.im_xcorr_coelution_score, 2.73205080756888) + // Test #3: IM Scores (Condition 2/2) dia_extract_window_ = 0.1; - IonMobilityScoring::driftScoring(drift_spectrum, transitions, scores, - drift_lower, drift_upper, drift_target, + IonMobilityScoring::driftScoring(sptrArr2, transitions, scores, + drift_target, im_range_1, dia_extract_window_, dia_extraction_ppm_, false, im_drift_extra_pcnt_); @@ -260,13 +271,13 @@ START_SECTION(([EXTRA] TEST_REAL_SIMILAR(scores.im_xcorr_shape_score, 0.833513903989399) TEST_REAL_SIMILAR(scores.im_xcorr_coelution_score, 0.910683602522959) - // deal with exactly one entry in mobilogram + // Test #4: deal with exactly one entry in mobilogram dia_extract_window_ = 0.3; - drift_lower = 1.0; - drift_upper = 1.1; drift_target = 1.05; - IonMobilityScoring::driftScoring(drift_spectrum, transitions, scores, - drift_lower, drift_upper, drift_target, + OpenMS::RangeMobility im_range_2(drift_target); + im_range_2.minSpanIfSingular(0.1); + IonMobilityScoring::driftScoring(sptrArr2, transitions, scores, + drift_target, im_range_2, dia_extract_window_, dia_extraction_ppm_, false, im_drift_extra_pcnt_); @@ -277,13 +288,14 @@ START_SECTION(([EXTRA] TEST_REAL_SIMILAR(scores.im_xcorr_shape_score, 0.0) // higher is better TEST_REAL_SIMILAR(scores.im_xcorr_coelution_score, 1.) // lower is better - // deal with one zero transitions + // Test #5: deal with one zero transitions dia_extract_window_ = 0.3; - drift_lower = 1.0; - drift_upper = 1.3; + OpenMS::RangeMobility im_range_3; + im_range_3.setMin(1.0); + im_range_3.setMax(1.3); drift_target = 1.1; - IonMobilityScoring::driftScoring(drift_spectrum, transitions, scores, - drift_lower, drift_upper, drift_target, + IonMobilityScoring::driftScoring(sptrArr2, transitions, scores, + drift_target, im_range_3, dia_extract_window_, dia_extraction_ppm_, false, im_drift_extra_pcnt_); @@ -294,12 +306,13 @@ START_SECTION(([EXTRA] TEST_REAL_SIMILAR(scores.im_xcorr_shape_score, 1.0/3) TEST_REAL_SIMILAR(scores.im_xcorr_coelution_score, 3.73205080756888) - // deal with all-zero transitions - drift_lower = 2.5; - drift_upper = 3.5; - drift_target = 3.0; - IonMobilityScoring::driftScoring(drift_spectrum, transitions, scores, - drift_lower, drift_upper, drift_target, + // Test #6: deal with all-zero transitions + // IM range from 2.5 to 3.5 + OpenMS::RangeMobility im_range_4(3.0); + im_range_4.minSpanIfSingular(1.); + + IonMobilityScoring::driftScoring(sptrArr2, transitions, scores, + drift_target, im_range_4, dia_extract_window_, dia_extraction_ppm_, false, im_drift_extra_pcnt_); @@ -313,7 +326,7 @@ START_SECTION(([EXTRA] END_SECTION START_SECTION([EXTRA] - static void driftScoringMS1(OpenSwath::SpectrumPtr spectrum, + static void driftScoringMS1(const SpectrumSequence& spectrum&, const std::vector & transitions, OpenSwath_Scores & scores, const double drift_lower, @@ -326,9 +339,10 @@ START_SECTION([EXTRA] { OpenSwath_Scores scores; - double drift_lower = 0.5; - double drift_upper = 1.5; + // IM range from 0.5 to 1.5 double drift_target = 1.0; + OpenMS::RangeMobility im_range(drift_target); + im_range.minSpanIfSingular(1.); double im_drift_extra_pcnt_ = 0.25; double dia_extract_window_ = 0.3; @@ -337,8 +351,11 @@ START_SECTION([EXTRA] OpenSwath::BinaryDataArrayPtr ion_mobility(new OpenSwath::BinaryDataArray); drift_spectrum->getDataArrays().push_back( ion_mobility ); - IonMobilityScoring::driftScoringMS1(drift_spectrum, transitions, scores, - drift_lower, drift_upper, drift_target, + std::vector sptrArr; + sptrArr.push_back(drift_spectrum); + + IonMobilityScoring::driftScoringMS1(sptrArr, transitions, scores, + drift_target, im_range, dia_extract_window_, dia_extraction_ppm_, false, im_drift_extra_pcnt_); @@ -347,19 +364,25 @@ START_SECTION([EXTRA] drift_spectrum->setMZArray( mass); drift_spectrum->setIntensityArray( intensity); - IonMobilityScoring::driftScoringMS1(drift_spectrum, transitions, scores, - drift_lower, drift_upper, drift_target, + std::vector sptrArr2; + sptrArr2.push_back(drift_spectrum); + + IonMobilityScoring::driftScoringMS1(sptrArr2, transitions, scores, + drift_target, im_range, dia_extract_window_, dia_extraction_ppm_, false, im_drift_extra_pcnt_); drift_spectrum = ms1spec; + std::vector sptrArr3; + sptrArr3.push_back(drift_spectrum); + TEST_EQUAL(drift_spectrum->getMZArray()->data.size(), 8) TEST_EQUAL(drift_spectrum->getMZArray()->data.size(), drift_spectrum->getIntensityArray()->data.size()) TEST_EQUAL(drift_spectrum->getMZArray()->data.size(), drift_spectrum->getDriftTimeArray()->data.size()) - IonMobilityScoring::driftScoringMS1(drift_spectrum, transitions, scores, - drift_lower, drift_upper, drift_target, + IonMobilityScoring::driftScoringMS1(sptrArr3, transitions, scores, + drift_target, im_range, dia_extract_window_, dia_extraction_ppm_, false, im_drift_extra_pcnt_); @@ -368,7 +391,7 @@ START_SECTION([EXTRA] END_SECTION START_SECTION(([EXTRA] - static void driftScoringMS1Contrast(OpenSwath::SpectrumPtr spectrum, OpenSwath::SpectrumPtr ms1spectrum, + static void driftScoringMS1Contrast(std::vector spectrum, OpenSwath::SpectrumPtr ms1spectrum, const std::vector & transitions, OpenSwath_Scores & scores, const double drift_lower, @@ -381,8 +404,9 @@ START_SECTION(([EXTRA] { OpenSwath_Scores scores; - double drift_lower = 0.5; - double drift_upper = 1.5; + // IM from 0.5 to 1.5 + OpenMS::RangeMobility im_range_1(1); + im_range_1.minSpanIfSingular(1.); double im_drift_extra_pcnt_ = 0.25; double dia_extract_window_ = 0.3; @@ -395,8 +419,14 @@ START_SECTION(([EXTRA] drift_spectrum_ms1->getDataArrays().push_back( ion_mobility ); drift_spectrum->getDataArrays().push_back( ion_mobility ); - IonMobilityScoring::driftScoringMS1Contrast(drift_spectrum, drift_spectrum_ms1, transitions, scores, - drift_lower, drift_upper, + std::vector sptrArr; + std::vector sptrArrMS1; + + sptrArr.push_back(drift_spectrum); + sptrArrMS1.push_back(drift_spectrum_ms1); + + IonMobilityScoring::driftScoringMS1Contrast(sptrArr, sptrArrMS1, transitions, scores, + im_range_1, dia_extract_window_, dia_extraction_ppm_, im_drift_extra_pcnt_); @@ -407,20 +437,34 @@ START_SECTION(([EXTRA] drift_spectrum_ms1->setMZArray( mass); drift_spectrum_ms1->setIntensityArray( intensity); - IonMobilityScoring::driftScoringMS1Contrast(drift_spectrum, drift_spectrum_ms1, transitions, scores, - drift_lower, drift_upper, + + std::vector sptrArr_2; + std::vector sptrArrMS1_2; + + sptrArr_2.push_back(drift_spectrum); + sptrArrMS1_2.push_back(drift_spectrum_ms1); + + IonMobilityScoring::driftScoringMS1Contrast(sptrArr_2, sptrArrMS1, transitions, scores, + im_range_1, dia_extract_window_, dia_extraction_ppm_, im_drift_extra_pcnt_); + + std::vector sptrArr_3; + std::vector sptrArrMS1_3; + drift_spectrum = spec; drift_spectrum_ms1 = ms1spec; + sptrArr_3.push_back(drift_spectrum); + sptrArrMS1_3.push_back(drift_spectrum_ms1); + TEST_EQUAL(drift_spectrum->getMZArray()->data.size(), 24) TEST_EQUAL(drift_spectrum->getMZArray()->data.size(), drift_spectrum->getIntensityArray()->data.size()) TEST_EQUAL(drift_spectrum->getMZArray()->data.size(), drift_spectrum->getDriftTimeArray()->data.size()) - - IonMobilityScoring::driftScoringMS1Contrast(drift_spectrum, drift_spectrum_ms1, transitions, scores, - drift_lower, drift_upper, + + IonMobilityScoring::driftScoringMS1Contrast(sptrArr_3, sptrArrMS1_3, transitions, scores, + im_range_1, dia_extract_window_, dia_extraction_ppm_, im_drift_extra_pcnt_); @@ -430,8 +474,8 @@ START_SECTION(([EXTRA] TEST_REAL_SIMILAR(scores.im_ms1_sum_contrast_shape, 0.56486260935015) dia_extract_window_ = 0.1; - IonMobilityScoring::driftScoringMS1Contrast(drift_spectrum, drift_spectrum_ms1, transitions, scores, - drift_lower, drift_upper, + IonMobilityScoring::driftScoringMS1Contrast(sptrArr_3, sptrArrMS1_3, transitions, scores, + im_range_1, dia_extract_window_, dia_extraction_ppm_, im_drift_extra_pcnt_); @@ -442,10 +486,13 @@ START_SECTION(([EXTRA] // deal with exactly one entry in mobilogram dia_extract_window_ = 0.3; - drift_lower = 1.0; - drift_upper = 1.1; - IonMobilityScoring::driftScoringMS1Contrast(drift_spectrum, drift_spectrum_ms1, transitions, scores, - drift_lower, drift_upper, + + // IM Span from 1.0 to 1.1 + OpenMS::RangeMobility im_range_2(1.05); + im_range_2.minSpanIfSingular(0.1); + + IonMobilityScoring::driftScoringMS1Contrast(sptrArr_3, sptrArrMS1_3, transitions, scores, + im_range_2, dia_extract_window_, dia_extraction_ppm_, im_drift_extra_pcnt_); @@ -456,10 +503,11 @@ START_SECTION(([EXTRA] // deal with one zero transitions dia_extract_window_ = 0.3; - drift_lower = 1.0; - drift_upper = 1.3; - IonMobilityScoring::driftScoringMS1Contrast(drift_spectrum, drift_spectrum_ms1, transitions, scores, - drift_lower, drift_upper, + //Im Span from 1.0 to 1.3 + OpenMS::RangeMobility im_range_3(1.15); + im_range_3.minSpanIfSingular(0.3); + IonMobilityScoring::driftScoringMS1Contrast(sptrArr_3, sptrArrMS1_3, transitions, scores, + im_range_3, dia_extract_window_, dia_extraction_ppm_, im_drift_extra_pcnt_); @@ -469,11 +517,12 @@ START_SECTION(([EXTRA] TEST_REAL_SIMILAR(scores.im_ms1_sum_contrast_shape, 0) // deal with all-zero transitions - drift_lower = 2.5; - drift_upper = 3.5; - - IonMobilityScoring::driftScoringMS1Contrast(drift_spectrum, drift_spectrum_ms1, transitions, scores, - drift_lower, drift_upper, + // IM span from 2.5 to 3.5 + OpenMS::RangeMobility im_range_4(3.); + im_range_4.minSpanIfSingular(1.); + + IonMobilityScoring::driftScoringMS1Contrast(sptrArr_3, sptrArrMS1_3, transitions, scores, + im_range_4, dia_extract_window_, dia_extraction_ppm_, im_drift_extra_pcnt_); diff --git a/src/tests/class_tests/openms/source/MRMFeatureScoring_test.cpp b/src/tests/class_tests/openms/source/MRMFeatureScoring_test.cpp index 98488282122..0c180a74754 100644 --- a/src/tests/class_tests/openms/source/MRMFeatureScoring_test.cpp +++ b/src/tests/class_tests/openms/source/MRMFeatureScoring_test.cpp @@ -123,6 +123,7 @@ START_SECTION((virtual void test_dia_scores())) MRMFeature mrmfeature = OpenSWATH_Test::createMockFeature(); int by_charge_state = 1; + RangeMobility empty_im_range; // find spectrum that is closest to the apex of the peak (set to 3120) using binary search MSSpectrum OpenMSspectrum = (*swath_map.RTBegin( 3120 )); @@ -171,7 +172,11 @@ START_SECTION((virtual void test_dia_scores())) // We have to reorder the transitions to make the tests work std::vector transitions = transition_group.getTransitions(); double isotope_corr = 0, isotope_overlap = 0; - diascoring.dia_isotope_scores(transitions, sptr, imrmfeature, isotope_corr, isotope_overlap); + + std::vector sptrArr; + sptrArr.push_back(sptr); + + diascoring.dia_isotope_scores(transitions, sptrArr, imrmfeature, empty_im_range, isotope_corr, isotope_overlap); delete imrmfeature; @@ -179,13 +184,13 @@ START_SECTION((virtual void test_dia_scores())) double ppm_score = 0, ppm_score_weighted = 0; std::vector ppm_errors; diascoring.dia_massdiff_score(transition_group.getTransitions(), - sptr, normalized_library_intensity, ppm_score, ppm_score_weighted, ppm_errors); + sptrArr, normalized_library_intensity, empty_im_range, ppm_score, ppm_score_weighted, ppm_errors); // Presence of b/y series score double bseries_score = 0, yseries_score = 0; String sequence = "SYVAWDR"; OpenMS::AASequence aas = AASequence::fromString(sequence); - diascoring.dia_by_ion_score(sptr, aas, by_charge_state, bseries_score, yseries_score); + diascoring.dia_by_ion_score(sptrArr, aas, by_charge_state, empty_im_range, bseries_score, yseries_score); TEST_REAL_SIMILAR(isotope_corr, 0.2866618 * transition_group.getTransitions().size() ) TEST_REAL_SIMILAR(isotope_corr, 0.85998565339479) @@ -206,7 +211,7 @@ START_SECTION((virtual void test_dia_scores())) // b/y series score with modifications bseries_score = 0, yseries_score = 0; aas.setModification(1, "Phospho" ); // modify the Y - diascoring.dia_by_ion_score(sptr, aas, by_charge_state, bseries_score, yseries_score); + diascoring.dia_by_ion_score(sptrArr, aas, by_charge_state, empty_im_range, bseries_score, yseries_score); TEST_EQUAL(bseries_score, 0) TEST_EQUAL(yseries_score, 1) } diff --git a/src/tests/class_tests/openms/source/OpenSwathScores_test.cpp b/src/tests/class_tests/openms/source/OpenSwathScores_test.cpp index 084ec728f38..dc0bcc1480d 100644 --- a/src/tests/class_tests/openms/source/OpenSwathScores_test.cpp +++ b/src/tests/class_tests/openms/source/OpenSwathScores_test.cpp @@ -51,7 +51,7 @@ START_SECTION((double get_quick_lda_score(double library_corr_, double library_n { OpenSwath_Scores scores; - TEST_REAL_SIMILAR( scores.get_quick_lda_score(1.0, 1.0, 1.0, 1.0, 1.0, 1.0), + TEST_REAL_SIMILAR( scores.get_quick_lda_score(1.0, 1.0, 1.0, 1.0, 1.0, 1.0), -0.5319046 + 2.1643962 + 8.0353047 + @@ -74,14 +74,14 @@ START_SECTION((double calculate_lda_prescore(const OpenSwath_Scores& scores) con scores.log_sn_score = 1.0; scores.elution_model_fit_score = 1.0; - TEST_REAL_SIMILAR( scores.calculate_lda_prescore(scores), - -0.34664267 + - 2.98700722 + - 7.05496384 + - 0.09445371 + - -5.71823862 + - -0.72989582 + - 1.88443209) + TEST_REAL_SIMILAR( scores.calculate_lda_prescore(scores), + -0.34664267 + + 2.98700722 + + 7.05496384 + + 0.09445371 + + -5.71823862 + + -0.72989582 + + 1.88443209) } @@ -100,13 +100,13 @@ START_SECTION((double calculate_lda_single_transition(const OpenSwath_Scores& sc scores.log_sn_score = 1.0; scores.elution_model_fit_score = 1.0; - TEST_REAL_SIMILAR( scores.calculate_lda_single_transition(scores), + TEST_REAL_SIMILAR( scores.calculate_lda_single_transition(scores), 7.05496384 + -0.72989582 + -1.08443209) } END_SECTION - + START_SECTION((double calculate_swath_lda_prescore(const OpenSwath_Scores& scores) const)) { OpenSwath_Scores scores; @@ -124,17 +124,17 @@ START_SECTION((double calculate_swath_lda_prescore(const OpenSwath_Scores& score scores.log_sn_score = 1.0; - TEST_REAL_SIMILAR( scores.calculate_swath_lda_prescore(scores), - -0.19011762 + - 2.47298914 + - 5.63906731 + - -0.62640133 + - 0.36006925 + - 0.08814003 + - 0.13978311 + - -1.16475032 + - -0.19267813 + - -0.61712054) + TEST_REAL_SIMILAR( scores.calculate_swath_lda_prescore(scores), + -0.19011762 + + 2.47298914 + + 5.63906731 + + -0.62640133 + + 0.36006925 + + 0.08814003 + + 0.13978311 + + -1.16475032 + + -0.19267813 + + -0.61712054) } END_SECTION diff --git a/src/tests/class_tests/openms/source/OpenSwathScoring_test.cpp b/src/tests/class_tests/openms/source/OpenSwathScoring_test.cpp index 0bfe1232068..b4a3182f17c 100644 --- a/src/tests/class_tests/openms/source/OpenSwathScoring_test.cpp +++ b/src/tests/class_tests/openms/source/OpenSwathScoring_test.cpp @@ -27,6 +27,28 @@ ThisShouldFailAtCompileTime = 0 using namespace OpenMS; using namespace std; +MSSpectrum generateImSpec(int k_min, int k_max, double rt) +{ + MSSpectrum imSpec; + DataArrays::FloatDataArray fda; + imSpec.setRT(rt); + + for (int k = k_min; k < k_max; k++) + { + Peak1D p; + p.setMZ(100. + k); + p.setIntensity(1.); + imSpec.push_back(p); + fda.push_back((double) k); + fda.setName("Ion Mobility"); + } + imSpec.getFloatDataArrays().push_back(fda); + + return imSpec; +} + + + START_TEST(OpenSwathScoring, "$Id$") ///////////////////////////////////////////////////////////// @@ -48,12 +70,12 @@ START_SECTION(~OpenSwathScoring()) } END_SECTION -START_SECTION(( void initialize(double rt_normalization_factor, int add_up_spectra, double spacing_for_spectra_resampling, const OpenSwath_Scores_Usage & su, const std::string& spectrum_addition_method) )) +START_SECTION(( void initialize(double rt_normalization_factor, int add_up_spectra, double spacing_for_spectra_resampling, const OpenSwath_Scores_Usage & su, const std::string& spectrum_addition_method, bool use_ms1_ion_mobility) )) { ptr = new OpenSwathScoring(); OpenSwath_Scores_Usage su; TEST_NOT_EQUAL(ptr, nullPointer) - ptr->initialize(100.0, 1, 0.01, 0.0, su, "simple"); + ptr->initialize(100.0, 1, 0.01, 0.0, su, "simple", true); delete ptr; } END_SECTION @@ -94,19 +116,22 @@ START_SECTION((void getNormalized_library_intensities_(const std::vector swath_maps, - double RT, int nr_spectra_to_add, double drift_lower, drift_upper))) + double RT, int nr_spectra_to_add, RangeMobility im_range))) { + + OpenMS::RangeMobility im_range_empty; // use this empty im range as input for all examples // test result for empty map { boost::shared_ptr swath_map (new PeakMap); OpenSwath::SpectrumAccessPtr swath_ptr = SimpleOpenMSSpectraFactory::getSpectrumAccessOpenMSPtr(swath_map); OpenSwathScoring sc; - OpenSwath::SpectrumPtr sp = sc.fetchSpectrumSwath(swath_ptr, 20.0, 1, 0, 0); + std::vector sp = sc.fetchSpectrumSwath(swath_ptr, 20.0, 1, im_range_empty); - TEST_EQUAL(sp->getMZArray()->data.empty(), true); + TEST_EQUAL(sp.empty(), true); } + // test result for map with single spectrum { PeakMap* eptr = new PeakMap; @@ -123,23 +148,26 @@ START_SECTION((OpenSwath::SpectrumPtr OpenSwathScoring::fetchSpectrumSwath(std:: TEST_EQUAL(swath_ptr->getNrSpectra(), 1) OpenSwathScoring sc; OpenSwath_Scores_Usage su; - sc.initialize(1.0, 1, 0.005, 0.0, su, "resample"); - OpenSwath::SpectrumPtr sp = sc.fetchSpectrumSwath(swath_ptr, 20.0, 1, 0, 0); + sc.initialize(1.0, 1, 0.005, 0.0, su, "resample", true); + + std::vector sp = sc.fetchSpectrumSwath(swath_ptr, 20.0, 1, im_range_empty); - TEST_EQUAL(sp->getMZArray()->data.size(), 1); - TEST_EQUAL(sp->getIntensityArray()->data.size(), 1); + TEST_EQUAL(sp.size(), 1); + TEST_EQUAL(sp[0]->getMZArray()->data.size(), 1); + TEST_EQUAL(sp[0]->getIntensityArray()->data.size(), 1); - TEST_REAL_SIMILAR(sp->getMZArray()->data[0], 20.0); - TEST_REAL_SIMILAR(sp->getIntensityArray()->data[0], 200.0); + TEST_REAL_SIMILAR(sp[0]->getMZArray()->data[0], 20.0); + TEST_REAL_SIMILAR(sp[0]->getIntensityArray()->data[0], 200.0); - sc.initialize(1.0, 1, 0.005, 0.0, su, "simple"); - sp = sc.fetchSpectrumSwath(swath_ptr, 20.0, 1, 0, 0); + sc.initialize(1.0, 1, 0.005, 0.0, su, "simple", true); + sp = sc.fetchSpectrumSwath(swath_ptr, 20.0, 1, im_range_empty); - TEST_EQUAL(sp->getMZArray()->data.size(), 1); - TEST_EQUAL(sp->getIntensityArray()->data.size(), 1); + TEST_EQUAL(sp.size(), 1); + TEST_EQUAL(sp[0]->getMZArray()->data.size(), 1); + TEST_EQUAL(sp[0]->getIntensityArray()->data.size(), 1); - TEST_REAL_SIMILAR(sp->getMZArray()->data[0], 20.0); - TEST_REAL_SIMILAR(sp->getIntensityArray()->data[0], 200.0); + TEST_REAL_SIMILAR(sp[0]->getMZArray()->data[0], 20.0); + TEST_REAL_SIMILAR(sp[0]->getIntensityArray()->data[0], 200.0); // delete eptr; } @@ -163,27 +191,37 @@ START_SECTION((OpenSwath::SpectrumPtr OpenSwathScoring::fetchSpectrumSwath(std:: TEST_EQUAL(swath_ptr->getNrSpectra(), 3) OpenSwathScoring sc; OpenSwath_Scores_Usage su; - sc.initialize(1.0, 1, 0.005, 0.0, su, "resample"); - OpenSwath::SpectrumPtr sp = sc.fetchSpectrumSwath(swath_ptr, 20.0, 3, 0, 0); + sc.initialize(1.0, 1, 0.005, 0.0, su, "resample", true); + std::vector sp = sc.fetchSpectrumSwath(swath_ptr, 20.0, 3, im_range_empty); + + TEST_EQUAL(sp.size(), 1); + TEST_EQUAL(sp[0]->getMZArray()->data.size(), 1); + TEST_EQUAL(sp[0]->getIntensityArray()->data.size(), 1); + + TEST_REAL_SIMILAR(sp[0]->getMZArray()->data[0], 20.0); + TEST_REAL_SIMILAR(sp[0]->getIntensityArray()->data[0], 600.0); + + sc.initialize(1.0, 1, 0.005, 0.0, su, "simple", true); + sp = sc.fetchSpectrumSwath(swath_ptr, 20.0, 3, im_range_empty); + TEST_EQUAL(sp.size(), 3); + + TEST_EQUAL(sp[0]->getMZArray()->data.size(), 1); + TEST_EQUAL(sp[0]->getIntensityArray()->data.size(), 1); - TEST_EQUAL(sp->getMZArray()->data.size(), 1); - TEST_EQUAL(sp->getIntensityArray()->data.size(), 1); - TEST_REAL_SIMILAR(sp->getMZArray()->data[0], 20.0); - TEST_REAL_SIMILAR(sp->getIntensityArray()->data[0], 600.0); + TEST_EQUAL(sp[1]->getMZArray()->data.size(), 1); + TEST_EQUAL(sp[1]->getIntensityArray()->data.size(), 1); - sc.initialize(1.0, 1, 0.005, 0.0, su, "simple"); - sp = sc.fetchSpectrumSwath(swath_ptr, 20.0, 3, 0, 0); - TEST_EQUAL(sp->getMZArray()->data.size(), 3); - TEST_EQUAL(sp->getIntensityArray()->data.size(), 3); + TEST_EQUAL(sp[2]->getMZArray()->data.size(), 1); + TEST_EQUAL(sp[2]->getIntensityArray()->data.size(), 1); - TEST_REAL_SIMILAR(sp->getMZArray()->data[0], 20.0); - TEST_REAL_SIMILAR(sp->getIntensityArray()->data[0], 200.0); - TEST_REAL_SIMILAR(sp->getMZArray()->data[1], 20.0); - TEST_REAL_SIMILAR(sp->getIntensityArray()->data[1], 200.0); - TEST_REAL_SIMILAR(sp->getMZArray()->data[2], 20.0); - TEST_REAL_SIMILAR(sp->getIntensityArray()->data[2], 200.0); + TEST_REAL_SIMILAR(sp[0]->getMZArray()->data[0], 20.0); + TEST_REAL_SIMILAR(sp[0]->getIntensityArray()->data[0], 200.0); + TEST_REAL_SIMILAR(sp[1]->getMZArray()->data[0], 20.0); + TEST_REAL_SIMILAR(sp[1]->getIntensityArray()->data[0], 200.0); + TEST_REAL_SIMILAR(sp[2]->getMZArray()->data[0], 20.0); + TEST_REAL_SIMILAR(sp[2]->getIntensityArray()->data[0], 200.0); // delete eptr; } @@ -220,34 +258,300 @@ START_SECTION((OpenSwath::SpectrumPtr OpenSwathScoring::fetchSpectrumSwath(std:: TEST_EQUAL(swath_ptr->getNrSpectra(), 4) OpenSwathScoring sc; OpenSwath_Scores_Usage su; - sc.initialize(1.0, 1, 0.005, 0.0, su, "resample"); - OpenSwath::SpectrumPtr sp = sc.fetchSpectrumSwath(swath_ptr, 20.0, 3, 0, 0); - - TEST_EQUAL(sp->getMZArray()->data.size(), 3); - TEST_EQUAL(sp->getIntensityArray()->data.size(), 3); - - TEST_REAL_SIMILAR(sp->getMZArray()->data[0], 20.0); - TEST_REAL_SIMILAR(sp->getIntensityArray()->data[0], 360.0); - TEST_REAL_SIMILAR(sp->getMZArray()->data[1], 20.005); - TEST_REAL_SIMILAR(sp->getIntensityArray()->data[1], 40.0); - TEST_REAL_SIMILAR(sp->getMZArray()->data[2], 250.0); - TEST_REAL_SIMILAR(sp->getIntensityArray()->data[2], 300.0); - - sc.initialize(1.0, 1, 0.005, 0.0, su, "simple"); - sp = sc.fetchSpectrumSwath(swath_ptr, 20.0, 3, 0, 0); - TEST_EQUAL(sp->getMZArray()->data.size(), 3); - TEST_EQUAL(sp->getIntensityArray()->data.size(), 3); - TEST_REAL_SIMILAR(sp->getMZArray()->data[0], 20.0); - TEST_REAL_SIMILAR(sp->getIntensityArray()->data[0], 200.0); - TEST_REAL_SIMILAR(sp->getMZArray()->data[1], 20.001); - TEST_REAL_SIMILAR(sp->getIntensityArray()->data[1], 200.0); - TEST_REAL_SIMILAR(sp->getMZArray()->data[2], 250.0); - TEST_REAL_SIMILAR(sp->getIntensityArray()->data[2], 300.0); + sc.initialize(1.0, 1, 0.005, 0.0, su, "resample", true); + std::vector sp = sc.fetchSpectrumSwath(swath_ptr, 20.0, 3, im_range_empty); + + TEST_EQUAL(sp.size(), 1); + TEST_EQUAL(sp[0]->getMZArray()->data.size(), 3); + TEST_EQUAL(sp[0]->getIntensityArray()->data.size(), 3); + + TEST_REAL_SIMILAR(sp[0]->getMZArray()->data[0], 20.0); + TEST_REAL_SIMILAR(sp[0]->getIntensityArray()->data[0], 360.0); + TEST_REAL_SIMILAR(sp[0]->getMZArray()->data[1], 20.005); + TEST_REAL_SIMILAR(sp[0]->getIntensityArray()->data[1], 40.0); + TEST_REAL_SIMILAR(sp[0]->getMZArray()->data[2], 250.0); + TEST_REAL_SIMILAR(sp[0]->getIntensityArray()->data[2], 300.0); + + // in simple method all 3 spectra should be returned + sc.initialize(1.0, 1, 0.005, 0.0, su, "simple", true); + sp = sc.fetchSpectrumSwath(swath_ptr, 20.0, 3, im_range_empty); + TEST_EQUAL(sp.size(), 3); + TEST_EQUAL(sp[0]->getMZArray()->data.size(), 1); + TEST_EQUAL(sp[1]->getMZArray()->data.size(), 1); + TEST_EQUAL(sp[2]->getMZArray()->data.size(), 1); + TEST_EQUAL(sp[0]->getIntensityArray()->data.size(), 1); + TEST_EQUAL(sp[1]->getIntensityArray()->data.size(), 1); + TEST_EQUAL(sp[2]->getIntensityArray()->data.size(), 1); + TEST_REAL_SIMILAR(sp[0]->getMZArray()->data[0], 20.001); + TEST_REAL_SIMILAR(sp[0]->getIntensityArray()->data[0], 200.0); + TEST_REAL_SIMILAR(sp[1]->getMZArray()->data[0], 20.0); + TEST_REAL_SIMILAR(sp[1]->getIntensityArray()->data[0], 200.0); + TEST_REAL_SIMILAR(sp[2]->getMZArray()->data[0], 250.0); + TEST_REAL_SIMILAR(sp[2]->getIntensityArray()->data[0], 300.0); } } END_SECTION -///////////////////////////////////////////////////////////// -///////////////////////////////////////////////////////////// -END_TEST +START_SECTION((OpenSwath::SpectrumPtr OpenSwathScoring::fetchSpectrumSwath(std::vector swath_maps, + double RT, int nr_spectra_to_add, RangeMobility im_range))- extra) +{ + // im range from 2-4 + OpenMS::RangeMobility im_range(3); // use this empty im range as input for all examples + im_range.minSpanIfSingular(2); // + + + // test result for map with single spectrum, should filter by IM because resampling is set + { + PeakMap* eptr = new PeakMap; + eptr->addSpectrum(generateImSpec(1,6,20.0)); + boost::shared_ptr swath_map (eptr); + OpenSwath::SpectrumAccessPtr swath_ptr = SimpleOpenMSSpectraFactory::getSpectrumAccessOpenMSPtr(swath_map); + TEST_EQUAL(swath_ptr->getNrSpectra(), 1); + + OpenSwathScoring sc; + OpenSwath_Scores_Usage su; + + // test resample - IM filtering should occur + { + sc.initialize(1.0, 1, 0.005, 0.0, su, "resample", true); + + SpectrumSequence sp = sc.fetchSpectrumSwath(swath_ptr, 20.0, 1, im_range); + + TEST_EQUAL(sp.size(), 1); + TEST_EQUAL(sp[0]->getMZArray()->data.size(), 3); + TEST_EQUAL(sp[0]->getIntensityArray()->data.size(), 3); + TEST_EQUAL(sp[0]->getDriftTimeArray()->data.size(), 3); + + TEST_REAL_SIMILAR(sp[0]->getMZArray()->data[0], 102.); + TEST_REAL_SIMILAR(sp[0]->getMZArray()->data[1], 103.); + TEST_REAL_SIMILAR(sp[0]->getMZArray()->data[2], 104.); + + TEST_REAL_SIMILAR(sp[0]->getIntensityArray()->data[0], 1.); + TEST_REAL_SIMILAR(sp[0]->getIntensityArray()->data[1], 1.); + TEST_REAL_SIMILAR(sp[0]->getIntensityArray()->data[2], 1.); + + TEST_REAL_SIMILAR(sp[0]->getDriftTimeArray()->data[0], 2.); + TEST_REAL_SIMILAR(sp[0]->getDriftTimeArray()->data[1], 3.); + TEST_REAL_SIMILAR(sp[0]->getDriftTimeArray()->data[2], 4.); + } + // test simple, since downstream functions are IM aware no filtering needs to occur. + { + sc.initialize(1.0, 1, 0.005, 0.0, su, "simple", true); + SpectrumSequence sp = sc.fetchSpectrumSwath(swath_ptr, 20.0, 1, im_range); + + TEST_EQUAL(sp.size(), 1); + TEST_EQUAL(sp[0]->getMZArray()->data.size(), 5); + TEST_EQUAL(sp[0]->getIntensityArray()->data.size(), 5); + + TEST_REAL_SIMILAR(sp[0]->getMZArray()->data[0], 101.); + TEST_REAL_SIMILAR(sp[0]->getMZArray()->data[1], 102.); + TEST_REAL_SIMILAR(sp[0]->getMZArray()->data[2], 103.); + TEST_REAL_SIMILAR(sp[0]->getMZArray()->data[3], 104.); + TEST_REAL_SIMILAR(sp[0]->getMZArray()->data[4], 105.); + + TEST_REAL_SIMILAR(sp[0]->getIntensityArray()->data[0], 1.); + TEST_REAL_SIMILAR(sp[0]->getIntensityArray()->data[1], 1.); + TEST_REAL_SIMILAR(sp[0]->getIntensityArray()->data[2], 1.); + TEST_REAL_SIMILAR(sp[0]->getIntensityArray()->data[3], 1.); + TEST_REAL_SIMILAR(sp[0]->getIntensityArray()->data[4], 1.); + + TEST_REAL_SIMILAR(sp[0]->getDriftTimeArray()->data[0], 1.); + TEST_REAL_SIMILAR(sp[0]->getDriftTimeArray()->data[1], 2.); + TEST_REAL_SIMILAR(sp[0]->getDriftTimeArray()->data[2], 3.); + TEST_REAL_SIMILAR(sp[0]->getDriftTimeArray()->data[3], 4.); + TEST_REAL_SIMILAR(sp[0]->getDriftTimeArray()->data[4], 5.); + } + // delete eptr; + } + + // Test result for 3 spectra + { + PeakMap* eptr = new PeakMap; + eptr->addSpectrum(generateImSpec(1,3,19.0)); + eptr->addSpectrum(generateImSpec(1,6,20.0)); + eptr->addSpectrum(generateImSpec(3,6,21.0)); + boost::shared_ptr swath_map (eptr); + OpenSwath::SpectrumAccessPtr swath_ptr = SimpleOpenMSSpectraFactory::getSpectrumAccessOpenMSPtr(swath_map); + TEST_EQUAL(swath_ptr->getNrSpectra(), 3); + + OpenSwathScoring sc; + OpenSwath_Scores_Usage su; + + // test resample - IM filtering should occur, also IM information is not needed so is cleared + { + sc.initialize(1.0, 1, 0.005, 0.0, su, "resample", true); + SpectrumSequence sp = sc.fetchSpectrumSwath(swath_ptr, 20.0, 3, im_range); + + TEST_EQUAL(sp.size(), 1); + + TEST_EQUAL(sp[0]->getMZArray()->data.size(), 3); + TEST_EQUAL(sp[0]->getIntensityArray()->data.size(), 3); + TEST_TRUE( (sp[0]->getDriftTimeArray() == nullptr ) ); // for resampling we do not use IM array + + TEST_REAL_SIMILAR(sp[0]->getMZArray()->data[0], 102.); + TEST_REAL_SIMILAR(sp[0]->getMZArray()->data[1], 103.); + TEST_REAL_SIMILAR(sp[0]->getMZArray()->data[2], 104.); + + TEST_REAL_SIMILAR(sp[0]->getIntensityArray()->data[0], 2.); + TEST_REAL_SIMILAR(sp[0]->getIntensityArray()->data[1], 2.); + TEST_REAL_SIMILAR(sp[0]->getIntensityArray()->data[2], 2.); + } + // test simple, since downstream functions are IM aware no filtering needs to occur. Should just return all the original spectra + { + sc.initialize(1.0, 1, 0.005, 0.0, su, "simple", true); + SpectrumSequence sp = sc.fetchSpectrumSwath(swath_ptr, 20.0, 3, im_range); + + //test sizing + TEST_EQUAL(sp.size(), 3); + TEST_EQUAL(sp[0]->getMZArray()->data.size(), 5); + TEST_EQUAL(sp[0]->getIntensityArray()->data.size(), 5); + TEST_EQUAL(sp[0]->getDriftTimeArray()->data.size(), 5); + + TEST_EQUAL(sp[1]->getMZArray()->data.size(), 2); + TEST_EQUAL(sp[1]->getIntensityArray()->data.size(), 2); + TEST_EQUAL(sp[1]->getDriftTimeArray()->data.size(), 2); + + TEST_EQUAL(sp[2]->getMZArray()->data.size(), 3); + TEST_EQUAL(sp[2]->getIntensityArray()->data.size(), 3); + TEST_EQUAL(sp[2]->getDriftTimeArray()->data.size(), 3); + + // Spectrum #1 + TEST_REAL_SIMILAR(sp[0]->getMZArray()->data[0], 101.); + TEST_REAL_SIMILAR(sp[0]->getMZArray()->data[1], 102.); + TEST_REAL_SIMILAR(sp[0]->getMZArray()->data[2], 103.); + TEST_REAL_SIMILAR(sp[0]->getMZArray()->data[3], 104.); + TEST_REAL_SIMILAR(sp[0]->getMZArray()->data[4], 105.); + + TEST_REAL_SIMILAR(sp[0]->getIntensityArray()->data[0], 1.); + TEST_REAL_SIMILAR(sp[0]->getIntensityArray()->data[1], 1.); + TEST_REAL_SIMILAR(sp[0]->getIntensityArray()->data[2], 1.); + TEST_REAL_SIMILAR(sp[0]->getIntensityArray()->data[3], 1.); + TEST_REAL_SIMILAR(sp[0]->getIntensityArray()->data[4], 1.); + + TEST_REAL_SIMILAR(sp[0]->getDriftTimeArray()->data[0], 1.); + TEST_REAL_SIMILAR(sp[0]->getDriftTimeArray()->data[1], 2.); + TEST_REAL_SIMILAR(sp[0]->getDriftTimeArray()->data[2], 3.); + TEST_REAL_SIMILAR(sp[0]->getDriftTimeArray()->data[3], 4.); + TEST_REAL_SIMILAR(sp[0]->getDriftTimeArray()->data[4], 5.); + + // Spectrum #2 + TEST_REAL_SIMILAR(sp[1]->getMZArray()->data[0], 101.); + TEST_REAL_SIMILAR(sp[1]->getMZArray()->data[1], 102.); + + TEST_REAL_SIMILAR(sp[1]->getIntensityArray()->data[0], 1.); + TEST_REAL_SIMILAR(sp[1]->getIntensityArray()->data[1], 1.); + + TEST_REAL_SIMILAR(sp[1]->getDriftTimeArray()->data[0], 1.); + TEST_REAL_SIMILAR(sp[1]->getDriftTimeArray()->data[1], 2.); + + // Spectrum #3 + TEST_REAL_SIMILAR(sp[2]->getMZArray()->data[0], 103.); + TEST_REAL_SIMILAR(sp[2]->getMZArray()->data[1], 104.); + + TEST_REAL_SIMILAR(sp[2]->getIntensityArray()->data[0], 1.); + TEST_REAL_SIMILAR(sp[2]->getIntensityArray()->data[1], 1.); + + TEST_REAL_SIMILAR(sp[2]->getDriftTimeArray()->data[0], 3.); + TEST_REAL_SIMILAR(sp[2]->getDriftTimeArray()->data[1], 4.); + } + // delete eptr; + } + + // test result for map with 4 spectra (select 3) + { + PeakMap* eptr = new PeakMap; + eptr->addSpectrum(generateImSpec(1,3,19.0)); + eptr->addSpectrum(generateImSpec(1,6,20.0)); + eptr->addSpectrum(generateImSpec(3,6,21.0)); + eptr->addSpectrum(generateImSpec(1,6,250)); + boost::shared_ptr swath_map (eptr); + OpenSwath::SpectrumAccessPtr swath_ptr = SimpleOpenMSSpectraFactory::getSpectrumAccessOpenMSPtr(swath_map); + TEST_EQUAL(swath_ptr->getNrSpectra(), 4); + + OpenSwathScoring sc; + OpenSwath_Scores_Usage su; + + // Test resampling, IM filtering should occur and the 4th spectrum should not be selected + { + sc.initialize(1.0, 1, 0.005, 0.0, su, "resample", true); + + SpectrumSequence sp = sc.fetchSpectrumSwath(swath_ptr, 20.0, 3, im_range); + + TEST_EQUAL(sp.size(), 1); + TEST_EQUAL(sp[0]->getMZArray()->data.size(), 3); + TEST_EQUAL(sp[0]->getIntensityArray()->data.size(), 3); + TEST_TRUE( (sp[0]->getDriftTimeArray() == nullptr ) ); // for resampling we do not use IM array + + TEST_REAL_SIMILAR(sp[0]->getMZArray()->data[0], 102.); + TEST_REAL_SIMILAR(sp[0]->getMZArray()->data[1], 103.); + TEST_REAL_SIMILAR(sp[0]->getMZArray()->data[2], 104.); + + TEST_REAL_SIMILAR(sp[0]->getIntensityArray()->data[0], 2.); + TEST_REAL_SIMILAR(sp[0]->getIntensityArray()->data[1], 2.); + TEST_REAL_SIMILAR(sp[0]->getIntensityArray()->data[2], 2.); + } + + // test simple, since downstream functions are IM aware no filtering needs to occur. Should just return all the original spectra, but the 4th spectrum should not be selected + { + sc.initialize(1.0, 1, 0.005, 0.0, su, "simple", true); + SpectrumSequence sp = sc.fetchSpectrumSwath(swath_ptr, 20.0, 3, im_range); + + //test sizing + TEST_EQUAL(sp.size(), 3); + TEST_EQUAL(sp[0]->getMZArray()->data.size(), 5); + TEST_EQUAL(sp[0]->getIntensityArray()->data.size(), 5); + TEST_EQUAL(sp[0]->getDriftTimeArray()->data.size(), 5); + + TEST_EQUAL(sp[1]->getMZArray()->data.size(), 2); + TEST_EQUAL(sp[1]->getIntensityArray()->data.size(), 2); + TEST_EQUAL(sp[1]->getDriftTimeArray()->data.size(), 2); + + TEST_EQUAL(sp[2]->getMZArray()->data.size(), 3); + TEST_EQUAL(sp[2]->getIntensityArray()->data.size(), 3); + TEST_EQUAL(sp[2]->getDriftTimeArray()->data.size(), 3); + + // Spectrum #1 + TEST_REAL_SIMILAR(sp[0]->getMZArray()->data[0], 101.); + TEST_REAL_SIMILAR(sp[0]->getMZArray()->data[1], 102.); + TEST_REAL_SIMILAR(sp[0]->getMZArray()->data[2], 103.); + TEST_REAL_SIMILAR(sp[0]->getMZArray()->data[3], 104.); + TEST_REAL_SIMILAR(sp[0]->getMZArray()->data[4], 105.); + + TEST_REAL_SIMILAR(sp[0]->getIntensityArray()->data[0], 1.); + TEST_REAL_SIMILAR(sp[0]->getIntensityArray()->data[1], 1.); + TEST_REAL_SIMILAR(sp[0]->getIntensityArray()->data[2], 1.); + TEST_REAL_SIMILAR(sp[0]->getIntensityArray()->data[3], 1.); + TEST_REAL_SIMILAR(sp[0]->getIntensityArray()->data[4], 1.); + + TEST_REAL_SIMILAR(sp[0]->getDriftTimeArray()->data[0], 1.); + TEST_REAL_SIMILAR(sp[0]->getDriftTimeArray()->data[1], 2.); + TEST_REAL_SIMILAR(sp[0]->getDriftTimeArray()->data[2], 3.); + TEST_REAL_SIMILAR(sp[0]->getDriftTimeArray()->data[3], 4.); + TEST_REAL_SIMILAR(sp[0]->getDriftTimeArray()->data[4], 5.); + + // Spectrum #2 + TEST_REAL_SIMILAR(sp[1]->getMZArray()->data[0], 101.); + TEST_REAL_SIMILAR(sp[1]->getMZArray()->data[1], 102.); + + TEST_REAL_SIMILAR(sp[1]->getIntensityArray()->data[0], 1.); + TEST_REAL_SIMILAR(sp[1]->getIntensityArray()->data[1], 1.); + + TEST_REAL_SIMILAR(sp[1]->getDriftTimeArray()->data[0], 1.); + TEST_REAL_SIMILAR(sp[1]->getDriftTimeArray()->data[1], 2.); + + // Spectrum #3 + TEST_REAL_SIMILAR(sp[2]->getMZArray()->data[0], 103.); + TEST_REAL_SIMILAR(sp[2]->getMZArray()->data[1], 104.); + + TEST_REAL_SIMILAR(sp[2]->getIntensityArray()->data[0], 1.); + TEST_REAL_SIMILAR(sp[2]->getIntensityArray()->data[1], 1.); + + TEST_REAL_SIMILAR(sp[2]->getDriftTimeArray()->data[0], 3.); + TEST_REAL_SIMILAR(sp[2]->getDriftTimeArray()->data[1], 4.); + } + // delete eptr; + } +} +END_SECTION +END_TEST diff --git a/src/tests/class_tests/openms/source/SpectrumAddition_test.cpp b/src/tests/class_tests/openms/source/SpectrumAddition_test.cpp index 609f498e3e1..c5a593cffa6 100644 --- a/src/tests/class_tests/openms/source/SpectrumAddition_test.cpp +++ b/src/tests/class_tests/openms/source/SpectrumAddition_test.cpp @@ -1,6 +1,6 @@ // Copyright (c) 2002-2023, The OpenMS Team -- EKU Tuebingen, ETH Zurich, and FU Berlin // SPDX-License-Identifier: BSD-3-Clause -// +// // -------------------------------------------------------------------------- // $Maintainer: Hannes Roest $ // $Authors: Hannes Roest $ @@ -22,6 +22,118 @@ START_TEST(SpectrumAddition, "$Id$") ///////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////// +START_SECTION((void sortSpectrumByMZ(OpenSwath::Spectrum& spec)) - No IM) +{ + OpenSwath::SpectrumPtr spec(new OpenSwath::Spectrum()); + OpenSwath::BinaryDataArrayPtr mass(new OpenSwath::BinaryDataArray); + OpenSwath::BinaryDataArrayPtr intensity(new OpenSwath::BinaryDataArray); + + // Intensity Sorted + std::vector intensSorted = { + 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 + }; + + // Mass Sorted + std::vector massSorted = { + 100, 101.5, 101.9, 102.0, 102.1, 102.11, 102.2, 102.25, 102.3, 102.4, 102.45 + }; + + // Intensity Not Sorted + std::vector intensNotSorted = { + 11, 4, 3, 5, 6, 7, 8, 9, 1, 2, 10 + }; + + // Mass Not Sorted + std::vector massNotSorted = { + 102.45, 102.0, 101.9, 102.1, 102.11, 102.2, 102.25, 102.3, 100, 101.5, 102.4 + }; + + // IM Not Sorted + std::vector imNotSorted = { + 11, 4, 3, 5, 6, 7, 8, 9, 1, 2, 10 + }; + + mass->data=massNotSorted; + intensity->data=intensNotSorted; + + spec->setMZArray(mass); + spec->setIntensityArray(intensity); + SpectrumAddition::sortSpectrumByMZ(*spec); + + TEST_EQUAL(spec->getMZArray()->data.size(), massSorted.size()); + TEST_EQUAL(spec->getIntensityArray()->data.size(), intensSorted.size()); + + for (size_t i=0; igetMZArray()->data[i]); + TEST_REAL_SIMILAR(intensSorted[i], spec->getIntensityArray()->data[i]); + } +} +END_SECTION + +START_SECTION((void sortSpectrumByMZ(OpenSwath::Spectrum& spec)) - With IM) +{ + OpenSwath::SpectrumPtr specIM(new OpenSwath::Spectrum()); + OpenSwath::BinaryDataArrayPtr mass(new OpenSwath::BinaryDataArray); + OpenSwath::BinaryDataArrayPtr intensity(new OpenSwath::BinaryDataArray); + OpenSwath::BinaryDataArrayPtr im(new OpenSwath::BinaryDataArray); + + // Intensity Sorted + std::vector intensSorted = { + 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 + }; + + // Mass Sorted + std::vector massSorted = { + 100, 101.5, 101.9, 102.0, 102.1, 102.11, 102.2, 102.25, 102.3, 102.4, 102.45 + }; + + // IM Sorted + std::vector imSorted = { + 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 + }; + + // Intensity Not Sorted + std::vector intensNotSorted = { + 11, 4, 3, 5, 6, 7, 8, 9, 1, 2, 10 + }; + + // Mass Not Sorted + std::vector massNotSorted = { + 102.45, 102.0, 101.9, 102.1, 102.11, 102.2, 102.25, 102.3, 100, 101.5, 102.4 + }; + + // IM Not Sorted + std::vector imNotSorted = { + 11, 4, 3, 5, 6, 7, 8, 9, 1, 2, 10 + }; + + // Create non sorted IM spectrum + mass->data=massNotSorted; + intensity->data=intensNotSorted; + im->data=imNotSorted; + + specIM->setMZArray(mass); + specIM->setIntensityArray(intensity); + specIM->setDriftTimeArray(im); + + mass->data=massNotSorted; + intensity->data=intensNotSorted; + im->data=imNotSorted; + + SpectrumAddition::sortSpectrumByMZ(*specIM); + TEST_EQUAL(specIM->getMZArray()->data.size(), massSorted.size()); + TEST_EQUAL(specIM->getIntensityArray()->data.size(), intensSorted.size()); + TEST_EQUAL(specIM->getDriftTimeArray()->data.size(), imSorted.size()); + for (size_t i=0; igetMZArray()->data[i]); + TEST_REAL_SIMILAR(intensSorted[i], specIM->getIntensityArray()->data[i]); + TEST_REAL_SIMILAR(imSorted[i], specIM->getDriftTimeArray()->data[i]); + } +} +END_SECTION + START_SECTION((static OpenSwath::SpectrumPtr addUpSpectra(std::vector< OpenSwath::SpectrumPtr > all_spectra, double sampling_rate, double filter_zeros)) ) { OpenSwath::SpectrumPtr spec1(new OpenSwath::Spectrum()); diff --git a/src/tests/class_tests/openms/source/SpectrumHelpers_test.cpp b/src/tests/class_tests/openms/source/SpectrumHelpers_test.cpp index 2aa42cffd9c..2eca9e850de 100644 --- a/src/tests/class_tests/openms/source/SpectrumHelpers_test.cpp +++ b/src/tests/class_tests/openms/source/SpectrumHelpers_test.cpp @@ -59,11 +59,18 @@ START_SECTION ( [EXTRA] testscorefunction) sptr->setMZArray( data1); sptr->setIntensityArray( data2); - double mzres, intensityres; - DIAHelpers::integrateWindow(sptr,499.,501.,mzres, intensityres); + std::vector sptrArr; + sptrArr.push_back(sptr); + + double mzres, intensityres, imres; + // mz range from 499 to 501 + RangeMZ mz_range(500.); + RangeMobility im_range_empty; + mz_range.minSpanIfSingular(2.); + DIAHelpers::integrateWindow(sptrArr, mzres, imres, intensityres, mz_range, im_range_empty); TEST_REAL_SIMILAR(mzres, 499.392014652015); - TEST_REAL_SIMILAR(intensityres,273 ); + TEST_REAL_SIMILAR(intensityres, 273 ); // >> exp = [240, 74, 39, 15, 0] > 121 / 500.338842975207 @@ -72,19 +79,19 @@ START_SECTION ( [EXTRA] testscorefunction) // >> pearsonr(exp, theo) // (0.99463189043051314, 0.00047175434098498532) // - DIAHelpers::integrateWindow(sptr,499.6,501.4,mzres, intensityres); + mz_range.setMin(499.6); + mz_range.setMax(501.4); + DIAHelpers::integrateWindow(sptrArr, mzres, imres, intensityres, mz_range, im_range_empty); - std::cout << "mz" << mzres << std::endl; - std::cout << "intensity" << intensityres << std::endl; TEST_REAL_SIMILAR(mzres, 500.338842975207); TEST_REAL_SIMILAR(intensityres,121 ); - std::vector wincenter, mzresv, intresv; + std::vector wincenter, mzresv, intresv, imresv; wincenter.push_back(300.); wincenter.push_back(200.); wincenter.push_back(500.); wincenter.push_back(600.); - DIAHelpers::integrateWindows(sptr,wincenter,0.5, intresv, mzresv); + DIAHelpers::integrateWindows(sptrArr,wincenter,0.5, intresv, mzresv, imresv, im_range_empty); TEST_REAL_SIMILAR(mzresv[0], 300); TEST_REAL_SIMILAR(intresv[0],0 ); TEST_REAL_SIMILAR(mzresv[1],200 ); diff --git a/src/tests/class_tests/openms/source/SwathMapMassCorrection_test.cpp b/src/tests/class_tests/openms/source/SwathMapMassCorrection_test.cpp index 807e73af7d9..debce777b99 100644 --- a/src/tests/class_tests/openms/source/SwathMapMassCorrection_test.cpp +++ b/src/tests/class_tests/openms/source/SwathMapMassCorrection_test.cpp @@ -1,6 +1,6 @@ // Copyright (c) 2002-2023, The OpenMS Team -- EKU Tuebingen, ETH Zurich, and FU Berlin // SPDX-License-Identifier: BSD-3-Clause -// +// // -------------------------------------------------------------------------- // $Maintainer: Hannes Roest $ // $Authors: Hannes Roest $ @@ -23,7 +23,7 @@ using namespace OpenMS; typedef OpenSwath::LightTransition TransitionType; -typedef std::map TransitionGroupMapPtrType; +typedef std::map TransitionGroupMapPtrType; OpenMS::MRMFeatureFinderScoring::TransitionGroupMapType getData() { @@ -211,7 +211,7 @@ START_SECTION( void correctMZ(OpenMS::MRMFeatureFinderScoring::TransitionGroupMa OpenSwath::SpectrumAccessPtr sptr2 = SimpleOpenMSSpectraFactory::getSpectrumAccessOpenMSPtr(exp2); OpenSwath::SwathMap map; - map.sptr = sptr; + map.sptr = sptr; map.lower = 400; map.upper = 425; map.center = 412.5; @@ -527,7 +527,7 @@ START_SECTION( void correctIM(const std::mapgetSpectrumById(0)->getMZArray()->data; - + TEST_REAL_SIMILAR(trafo_result.apply(10), 0.889721627408994) TEST_REAL_SIMILAR(trafo_result.apply(20), 10.0974304068522) @@ -619,7 +619,7 @@ START_SECTION( void correctIM(const std::mapgetSpectrumById(0)->getMZArray()->data; - + TEST_REAL_SIMILAR(trafo_result.apply(10), 0.889721627408994) TEST_REAL_SIMILAR(trafo_result.apply(20), 10.0974304068522) @@ -643,7 +643,7 @@ START_SECTION( void correctIM(const std::mapgetSpectrumById(0)->getMZArray()->data; - + // only got a single peptide, so regression is only intercept TEST_REAL_SIMILAR(trafo_result.apply(10), 11) TEST_REAL_SIMILAR(trafo_result.apply(20), 11) @@ -672,6 +672,7 @@ START_SECTION( void correctIM(const std::map + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eJxzUGGAAy4g7sABALAcDS8= + + + + + + eJxzONTFAAJbDk8/lXoh6+rVmyfuRj6c9WTG89hXH95Evp/zce/ng19Xf6/5afH7wZ/Sf7/+l/578Mfid83P1d8Pft37ec7HyPcf3sS+mvF81pPIhyfuXr2ZdTX1wvRTWw7b7nbcMHnB114EBABbE0tI + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eJxzUGEAAyYg5gHiDhwAALCUDTM= + + + + + + eJwBbACT/0DChoAAAAAAAACXidGisLFFvH7EN8vo0NbVMNoV3pvh0uTI54XqE+1277Tx0fPR9bf3hvk+++P8d/76/3f+4/w++4b5t/fR9dHztPF27xPtherI59Lkm+EV3jDa1tXo0DfLfsRFvLCx0aKXiW8HSHA= + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eJxzUGEAAyYg5gHiDhwAALCUDTM= + + + + + + eJwBbACT/0DA1oAAAAAAAAD9kx+rtLhcwtrJ+s8p1abZnN0n4VvkSef66Xjsy+738AHz7fS/9nn4Hfqu+y39nP78/5z+Lf2u+x36efi/9u30AfP38MvueOz66UnnW+Qn4Zzdptkp1frP2slcwrS4H6v9k4dPSJg= + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eJxzUGEAAyYg5gHiDhwAALCUDTM= + + + + + + eJwBbACT/0C/+AAAAAAAAABRmVivQ7xvxY7MYNJM15DbUt+w4rvlg+gS63DtpO+08aTzd/Uy99X4ZPrh+039qf73/6n+Tf3h+2T61fgy93f1pPO08aTvcO0S64Pou+Ww4lLfkNtM12DSjsxvxUO8WK9RmdU9RXI= + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eJxzUGEAAyYg5gHiDhwAALCUDTM= + + + + + + eJwBbACT/0C+3QAAAAAAAADUnCCynL54x1nO+NO52NfceeC546nmWOnQ6xruO/A48hf02/WG9xv5nPoM/Gv9u/79/7v+a/0M/Jz6G/mG99v1F/Q48jvwGu7Q61jpqea543ng19y52PjTWc54x5y+ILLUnOyaRgo= + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eJxzUGEAAyYg5gHiDhwAALCUDTM= + + + + + + eJwBbACT/0C+DgAAAAAAAABknya0T8DxyKTPHdW/2cHdS+F15FLn7+lX7JHupPCU8mb0Hva+90j5wPol/Hv9w/79/8P+e/0l/MD6SPm+9x72ZvSU8qTwke5X7O/pUud15Evhwd2/2R3VpM/xyE/AJrRkn51kRV0= + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eJxzUGGAAy4g7sABALAcDS8= + + + + + + eJxzONTFAAJbDk8/lXoh6+rVmyfuRj6c9WTG89hXH95Evp/zce/ng19Xf6/5afH7wZ/Sf7/+l/578Mfid83P1d8Pft37ec7HyPcf3sS+mvF81pPIhyfuXr2ZdTX1wvRTWw7b7nbcMHnB114EBABbE0tI + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eJxzUGGAAy4g7sABALAcDS8= + + + + + + eJxzONTFAAJbDk8/lXoh6+rVmyfuRj6c9WTG89hXH95Evp/zce/ng19Xf6/5afH7wZ/Sf7/+l/578Mfid83P1d8Pft37ec7HyPcf3sS+mvF81pPIhyfuXr2ZdTX1wvRTWw7b7nbcMHnB114EBABbE0tI + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eJxzUGEAA2Yg5gXiDhwAALDQDTU= + + + + + + eJwBbACT/0DBLAAAAAAAAADpkXipULcmwcrICc9R1OXY7tyL4NDjzOaL6RbsdO6r8L/ytfSQ9lP4APqY+x/9lf78/5X+H/2Y+wD6U/iQ9rX0v/Kr8HTuFuyL6czm0OOL4O7c5dhR1AnPysgmwVC3eKnpkbe7Rbs= + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eJxzUGEAA2Yg5gXiDhwAALDQDTU= + + + + + + eJwBbACT/0DChoAAAAAAAACXidGisLFFvH7EN8vo0NbVMNoV3pvh0uTI54XqE+1277Tx0fPR9bf3hvk+++P8d/76/3f+4/w++4b5t/fR9dHztPF27xPtherI59Lkm+EV3jDa1tXo0DfLfsRFvLCx0aKXiW8HSHA= + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eJxzUGEAA2Yg5gXiDhwAALDQDTU= + + + + + + eJwBbACT/0C/bAAAAAAAAAANm7iwa71xxnHNKdMA2DHc5N8z4zHm7Ohw68Tt7u/18d3zqPVb9/f4gPr2+1v9sv76/7L+W/32+4D69/hb96j13fP18e7vxO1w6+zoMeYz4+TfMdwA2CnTcc1xxmu9uLANmzHUR+k= + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eJxzUGEAA2Yg5gXiDhwAALDQDTU= + + + + + + eJwBbACT/0DA1oAAAAAAAAD9kx+rtLhcwtrJ+s8p1abZnN0n4VvkSef66Xjsy+738AHz7fS/9nn4Hfqu+y39nP78/5z+Lf2u+x36efi/9u30AfP38MvueOz66UnnW+Qn4Zzdptkp1frP2slcwrS4H6v9k4dPSJg= + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eJxzUGEAA2Yg5gXiDhwAALDQDTU= + + + + + + eJwBbACT/0C97wAAAAAAAADGn3O0kMApydXPSdXl2eTdauGR5GvnBupr7KLus/Ch8nL0KPbG90/5xfop/H79xP78/8T+fv0p/MX6T/nG9yj2cvSh8rPwou5r7Abqa+eR5Grh5N3l2UnV1c8pyZDAc7TGnz30R5g= + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eJxzUGEAA2Yg5gXiDhwAALDQDTU= + + + + + + eJwBbACT/0C/+AAAAAAAAABRmVivQ7xvxY7MYNJM15DbUt+w4rvlg+gS63DtpO+08aTzd/Uy99X4ZPrh+039qf73/6n+Tf3h+2T61fgy93f1pPO08aTvcO0S64Pou+Ww4lLfkNtM12DSjsxvxUO8WK9RmdU9RXI= + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eJxzUGEAA2Yg5gXiDhwAALDQDTU= + + + + + + eJwBbACT/0C89QAAAAAAAADbouC2mcLrymDRp9Yd2/reY+Jw5TLot+oI7S3vLfEL88z0dPYF+IH56vpD/I39yP73/8j+jf1D/Or6gfkF+HT2zPQL8y3xLe8I7bfqMuhw5WPi+t4d26fWYNHrypnC4Lbbov/sRlY= + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eJxzUGEAA2Yg5gXiDhwAALDQDTU= + + + + + + eJwBbACT/0C+3QAAAAAAAADUnCCynL54x1nO+NO52NfceeC546nmWOnQ6xruO/A48hf02/WG9xv5nPoM/Gv9u/79/7v+a/0M/Jz6G/mG99v1F/Q48jvwGu7Q61jpqea543ng19y52PjTWc54x5y+ILLUnOyaRgo= + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eJxzUGEAA2Yg5gXiDhwAALDQDTU= + + + + + + eJwBbACT/0C8PwAAAAAAAAAkpa24HsQ8zInSrtcJ3M7fIeMb5szoQeuD7ZvvjvFh8xf1tPY7+K75D/tf/KD91P77/9T+oP1f/A/7rvk7+LT2F/Vh847xm++D7UHrzOgb5iHjzt8J3K7XidI8zB7ErbgkpRY4QnY= + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eJxzUGEAA2Yg5gXiDhwAALDQDTU= + + + + + + eJwBbACT/0C+DgAAAAAAAABknya0T8DxyKTPHdW/2cHdS+F15FLn7+lX7JHupPCU8mb0Hva+90j5wPol/Hv9w/79/8P+e/0l/MD6SPm+9x72ZvSU8qTwke5X7O/pUud15Evhwd2/2R3VpM/xyE/AJrRkn51kRV0= + + + + + + + + + 3268 + 5564 + 8046 + 10528 + 13010 + 15492 + 17974 + 20270 + 22566 + 25051 + 27536 + 30022 + 32508 + 34994 + 37480 + 39966 + 42452 + 44938 + + +47458 +0 + \ No newline at end of file diff --git a/src/tests/topp/OpenSwathWorkflow_23_b_output.featureXML b/src/tests/topp/OpenSwathWorkflow_23_b_output.featureXML new file mode 100644 index 00000000000..f05838897ee --- /dev/null +++ b/src/tests/topp/OpenSwathWorkflow_23_b_output.featureXML @@ -0,0 +1,493 @@ + + + + + + + + + + 25.200000476837158 + 412.5 + 3.726e05 + 0.0 + 0.0 + 0.70849 + 0 + + + 25.200000476837158 + 100.0 + 2.484e04 + 0.0 + 0.0 + 0.0 + 0 + + + + + + + + + + 25.200000476837158 + 101.0 + 4.968e04 + 0.0 + 0.0 + 0.0 + 0 + + + + + + + + + + 25.200000476837158 + 102.0 + 7.452e04 + 0.0 + 0.0 + 0.0 + 0 + + + + + + + + + + 25.200000476837158 + 103.0 + 9.936e04 + 0.0 + 0.0 + 0.0 + 0 + + + + + + + + + + 25.200000476837158 + 104.0 + 1.242e05 + 0.0 + 0.0 + 0.0 + 0 + + + + + + + + + + 25.200000476837158 + 412.5 + 2.3806e04 + 0.0 + 0.0 + 0.0 + 2 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 25.299999523162843 + 417.5 + 6.42735e05 + 0.0 + 0.0 + 0.702851 + 0 + + + 25.299999523162843 + 100.010000000000005 + 4.2849e04 + 0.0 + 0.0 + 0.0 + 0 + + + + + + + + + + 25.299999523162843 + 101.010000000000005 + 8.5698e04 + 0.0 + 0.0 + 0.0 + 0 + + + + + + + + + + 25.299999523162843 + 102.010000000000005 + 1.28547e05 + 0.0 + 0.0 + 0.0 + 0 + + + + + + + + + + 25.299999523162843 + 103.010000000000005 + 1.71396e05 + 0.0 + 0.0 + 0.0 + 0 + + + + + + + + + + 25.299999523162843 + 104.010000000000005 + 2.14245e05 + 0.0 + 0.0 + 0.0 + 0 + + + + + + + + + + 25.299999523162843 + 417.5 + 2.4021e04 + 0.0 + 0.0 + 0.0 + 2 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 25.300000476837159 + 422.5 + 3.726e05 + 0.0 + 0.0 + 0.702851 + 0 + + + 25.300000476837159 + 100.019999999999996 + 2.484e04 + 0.0 + 0.0 + 0.0 + 0 + + + + + + + + + + 25.300000476837159 + 101.019999999999996 + 4.968e04 + 0.0 + 0.0 + 0.0 + 0 + + + + + + + + + + 25.300000476837159 + 102.019999999999996 + 7.452e04 + 0.0 + 0.0 + 0.0 + 0 + + + + + + + + + + 25.300000476837159 + 103.019999999999996 + 9.936e04 + 0.0 + 0.0 + 0.0 + 0 + + + + + + + + + + 25.300000476837159 + 104.019999999999996 + 1.242e05 + 0.0 + 0.0 + 0.0 + 0 + + + + + + + + + + 25.300000476837159 + 422.5 + 2.4021e04 + 0.0 + 0.0 + 0.0 + 2 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/topp/OpenSwathDIAPreScoring.cpp b/src/topp/OpenSwathDIAPreScoring.cpp index 704963abe65..cca63491578 100644 --- a/src/topp/OpenSwathDIAPreScoring.cpp +++ b/src/topp/OpenSwathDIAPreScoring.cpp @@ -141,7 +141,8 @@ class DIAPreScoring : swath_map); OpenSwath::IDataFrameWriter* dfw = new OpenSwath::CSVWriter(fname); OpenMS::DiaPrescore dp; - dp.operator()(spectrumAccess, transition_exp_used, dfw); + OpenMS::RangeMobility im_range; // create empty IM range object + dp.operator()(spectrumAccess, transition_exp_used, im_range, dfw); //note IM not supported here yet delete dfw; } //end of for loop return EXECUTION_OK; diff --git a/src/topp/OpenSwathWorkflow.cpp b/src/topp/OpenSwathWorkflow.cpp index 78518a9cc7d..6893ec89f91 100644 --- a/src/topp/OpenSwathWorkflow.cpp +++ b/src/topp/OpenSwathWorkflow.cpp @@ -752,9 +752,11 @@ class TOPPOpenSwathWorkflow ChromExtractParams cp_ms1 = cp; cp_ms1.mz_extraction_window = getDoubleOption_("mz_extraction_window_ms1"); cp_ms1.ppm = getStringOption_("mz_extraction_window_ms1_unit") == "ppm"; - cp_ms1.im_extraction_window = getDoubleOption_("im_extraction_window_ms1"); + cp_ms1.im_extraction_window = (use_ms1_im) ? getDoubleOption_("im_extraction_window_ms1") : -1; Param feature_finder_param = getParam_().copy("Scoring:", true); + feature_finder_param.setValue("use_ms1_ion_mobility", getStringOption_("use_ms1_ion_mobility")); + Param tsv_reader_param = getParam_().copy("Library:", true); if (use_emg_score) {