Repository: matchms/matchms
Pull request #663 (2024-05-27 12:43:52)

Remove noise below frequent intensities


Patch:
diff --git a/matchms/filtering/__init__.py b/matchms/filtering/__init__.py
index b7018091..82b92bb1 100644
--- a/matchms/filtering/__init__.py
+++ b/matchms/filtering/__init__.py
@@ -145,6 +145,8 @@ from matchms.filtering.peak_processing.normalize_intensities import \
     normalize_intensities
 from matchms.filtering.peak_processing.reduce_to_number_of_peaks import \
     reduce_to_number_of_peaks
+from matchms.filtering.peak_processing.remove_noise_below_frequent_intensities import \
+    remove_noise_below_frequent_intensities
 from matchms.filtering.peak_processing.remove_peaks_around_precursor_mz import \
     remove_peaks_around_precursor_mz
 from matchms.filtering.peak_processing.remove_peaks_outside_top_k import \
@@ -194,6 +196,7 @@ __all__ = [
     "repair_adduct_and_parent_mass_based_on_smiles",
     "repair_adduct_based_on_parent_mass",
     "require_correct_ionmode",
+    "remove_noise_below_frequent_intensities",
     "repair_inchi_inchikey_smiles",
     "repair_parent_mass_is_molar_mass",
     "repair_parent_mass_match_smiles_wrapper",
diff --git a/matchms/filtering/default_pipelines.py b/matchms/filtering/default_pipelines.py
index de255496..fd958579 100644
--- a/matchms/filtering/default_pipelines.py
+++ b/matchms/filtering/default_pipelines.py
@@ -49,6 +49,7 @@ CLEAN_PEAKS = [(msfilters.select_by_mz, {"mz_from": 0, "mz_to": 1000}),
                (msfilters.reduce_to_number_of_peaks, {"n_max": 1000}),
                (msfilters.require_minimum_number_of_high_peaks, {"no_peaks": 5, "intensity_percent": 2.0}),
                msfilters.remove_profiled_spectra,
+               msfilters.remove_noise_below_frequent_intensities
                ]
 # These filters are in None of the above pipelines
 OTHER_FILTERS = [matchms.filtering.metadata_processing.require_precursor_mz.require_precursor_below_mz,
diff --git a/matchms/filtering/filter_order.py b/matchms/filtering/filter_order.py
index 5d9a0a9b..a1fdc236 100644
--- a/matchms/filtering/filter_order.py
+++ b/matchms/filtering/filter_order.py
@@ -45,6 +45,7 @@ ALL_FILTERS = [msfilters.make_charge_int,
                msfilters.require_matching_adduct_and_ionmode,
                msfilters.remove_profiled_spectra,
                msfilters.normalize_intensities,
+               msfilters.remove_noise_below_frequent_intensities,
                msfilters.select_by_intensity,
                msfilters.select_by_mz,
                msfilters.select_by_relative_intensity,
diff --git a/matchms/filtering/peak_processing/remove_noise_below_frequent_intensities.py b/matchms/filtering/peak_processing/remove_noise_below_frequent_intensities.py
new file mode 100644
index 00000000..562793eb
--- /dev/null
+++ b/matchms/filtering/peak_processing/remove_noise_below_frequent_intensities.py
@@ -0,0 +1,55 @@
+import logging
+import numpy as np
+from matchms.Fragments import Fragments
+from matchms.Spectrum import Spectrum
+
+
+logger = logging.getLogger("matchms")
+
+
+def remove_noise_below_frequent_intensities(spectrum: Spectrum,
+                                            min_count_of_frequent_intensities: int = 5,
+                                            noise_level_multiplier: float = 2.0):
+    """Removes noise if intensities exactly match frequently
+
+    When no noise filtering has been applied to a spectrum, many spectra show repeating intensities.
+    From all intensities that repeat more than min_count_of_frequent_intensities the highest is selected.
+    The noise level is set to this intensity * noise_level_multiplier. All fragments with an intensity below the noise
+    level are removed.
+
+    This filter was suggested by Tytus Mak.
+
+    Parameters
+    ----------
+    spectrum
+        Input spectrum.
+    min_count_of_frequent_intensities:
+        Minimum number of repeating intensities.
+    noise_level_multiplier:
+        From all intensities that repeat more than min_count_of_frequent_intensities the highest is selected.
+    The noise level is set to this intensity * noise_level_multiplier.
+    """
+    if spectrum is None:
+        return None
+    spectrum = spectrum.clone()
+
+    highest_frequent_peak = _select_highest_frequent_peak(spectrum.intensities, min_count_of_frequent_intensities)
+    if highest_frequent_peak != -1:
+        noise_threshold = highest_frequent_peak * noise_level_multiplier
+        peaks_to_keep = spectrum.intensities > noise_threshold
+        new_mzs, new_intensities = spectrum.mz[peaks_to_keep], spectrum.intensities[peaks_to_keep]
+        spectrum.peaks = Fragments(mz=new_mzs, intensities=new_intensities)
+        logger.info("Fragments removed with intensity below %s", noise_threshold)
+    return spectrum
+
+
+def _select_highest_frequent_peak(intensities,
+                                  min_count_of_frequent_intensities=5):
+    unique_values, counts = np.unique(intensities, return_counts=True)
+    mask = counts >= min_count_of_frequent_intensities
+    filtered_values = unique_values[mask]
+
+    # Return the highest value from the filtered values, or -1 if no values meet the criteria
+    if filtered_values.size > 0:
+        return filtered_values.max()
+    return -1


Test patch:
diff --git a/tests/filtering/test_all_filter_order.py b/tests/filtering/test_all_filter_order.py
index 134364c7..264b2dfe 100644
--- a/tests/filtering/test_all_filter_order.py
+++ b/tests/filtering/test_all_filter_order.py
@@ -25,7 +25,7 @@ DERIVE_ANNOTATION_FILTERS = [msfilters.derive_smiles_from_inchi,
     [DERIVE_ANNOTATION_FILTERS + [msfilters.repair_inchi_inchikey_smiles, msfilters.clean_adduct,
                                   msfilters.derive_annotation_from_compound_name],
      REPAIR_PARENT_MASS_SMILES_FILTERS],
-    [[msfilters.add_precursor_mz,], [msfilters.require_precursor_mz,]],
+    [[msfilters.add_precursor_mz, ], [msfilters.require_precursor_mz, ]],
     # Since pubchem lookup checks if annotation is complete.
     # So deriving inchi and inchikey from smiles, should happen first.
     [[msfilters.repair_inchi_inchikey_smiles, msfilters.add_parent_mass],
@@ -39,6 +39,12 @@ DERIVE_ANNOTATION_FILTERS = [msfilters.derive_smiles_from_inchi,
     [[msfilters.derive_adduct_from_name, ], [msfilters.clean_adduct]],
     [[msfilters.derive_annotation_from_compound_name, ], DERIVE_ANNOTATION_FILTERS],
     [[msfilters.derive_formula_from_name, ], [msfilters.require_formula]],
+    [[msfilters.remove_profiled_spectra, ], [msfilters.remove_peaks_around_precursor_mz]],
+    [[msfilters.remove_noise_below_frequent_intensities],
+     [msfilters.select_by_intensity, msfilters.select_by_mz, msfilters.select_by_relative_intensity,
+      msfilters.remove_peaks_around_precursor_mz, msfilters.remove_peaks_outside_top_k,
+      msfilters.reduce_to_number_of_peaks, msfilters.require_minimum_number_of_peaks,
+      msfilters.require_minimum_number_of_high_peaks]],
     [[msfilters.remove_profiled_spectra,], [msfilters.remove_peaks_around_precursor_mz]],
     [[msfilters.clean_adduct, msfilters.derive_adduct_from_name, msfilters.repair_adduct_based_on_parent_mass,
       msfilters.repair_adduct_and_parent_mass_based_on_smiles, msfilters.add_precursor_mz,
@@ -68,6 +74,7 @@ def test_all_filters_is_complete():
 
     This is important, since performing tests in the wrong order can make some filters useless.
     """
+
     def get_functions_from_file(file_path):
         """Gets all python functions in a file"""
         with open(file_path, 'r', encoding="utf-8") as file:
diff --git a/tests/filtering/test_remove_noise_below_frequent_intensities.py b/tests/filtering/test_remove_noise_below_frequent_intensities.py
new file mode 100644
index 00000000..fc2fbd10
--- /dev/null
+++ b/tests/filtering/test_remove_noise_below_frequent_intensities.py
@@ -0,0 +1,25 @@
+import numpy as np
+import pytest
+from matchms.filtering import remove_noise_below_frequent_intensities
+from matchms.Spectrum import Spectrum
+
+
+@pytest.mark.parametrize("mz, intensities, min_count_of_frequent_intensities, noise_level_multiplier, "
+                         "expected_number_of_peaks", [
+                             [np.array([10, 20, 30, 40, 50, 60, 70, 80, 90], dtype="float"),
+                              np.array([20, 20, 20, 20, 20, 500, 10, 5, 20.1], dtype="float"), 5, 2.0, 1],
+                             [np.array([10, 20, 30, 40, 50, 60, 70, 80, 90], dtype="float"),
+                              np.array([20, 20, 20, 20, 20, 500, 10, 5, 100], dtype="float"), 5, 5.0, 1], # test noise level multiplier
+                             [np.array([10, 20, 30, 40, 50, 60, 70, 80, 90], dtype="float"),
+                              np.array([20, 20, 20, 20, 20, 500, 10, 5, 100], dtype="float"), 10, 5.0, 9], # test expected_number_of_peaks
+                         ])
+def test_remove_noise_below_frequent_intensities(mz, intensities, min_count_of_frequent_intensities,
+                                                 noise_level_multiplier,
+                                                 expected_number_of_peaks):
+    spectrum_in = Spectrum(mz, intensities)
+    spectrum_out = remove_noise_below_frequent_intensities(
+        spectrum_in,
+        min_count_of_frequent_intensities=min_count_of_frequent_intensities,
+        noise_level_multiplier=noise_level_multiplier)
+
+    assert len(spectrum_out.peaks.mz) == expected_number_of_peaks