Repository: matchms/matchms
Pull request #37 (2020-05-20 15:15:21)

236 modified cosine score
Originally here https://github.com/matchms/matchms-backup/pull/237


 -  Added ``ModifiedCosine`` class to similarities.
 -  Added unit test.
 -  Added ``numba`` to environment.


cwmeijer: This issue is marked for needing discussion but it has a PR already... Simple comment anyway:
Yay type checking! I never used it but have been planning to. User facing API is a great place to start!
jspaaks: The ``needs_discussion`` label was added March 23 though.
LourensVeen: When you start using pytest-mypy, you'll probably see some error messages related to libraries that don't have type annotations. See https://github.com/yatiml/yatiml/blob/develop/setup.cfg for an example of how to configure it to ignore those libraries. I'm using `pytest-mypy>=0.4.0` there to avoid a compatibility issue with recent versions of mypy, you may want to do the same. Also, don't forget to add a ``py.typed`` file (see https://www.python.org/dev/peps/pep-0561/).

I'm happy to help if you get stuck somewhere, just let me know.
jspaaks: PR #154

Patch:
diff --git a/conda/environment-dev.yml b/conda/environment-dev.yml
index 6c4d8d9a..eb4bdb45 100644
--- a/conda/environment-dev.yml
+++ b/conda/environment-dev.yml
@@ -6,6 +6,7 @@ channels:
 dependencies:
   - gensim >=3.8.0
   - matplotlib
+  - numba >=0.47
   - numpy
   - pip
   - pyteomics >=4.2
diff --git a/conda/environment.yml b/conda/environment.yml
index c4c51392..5af9b494 100644
--- a/conda/environment.yml
+++ b/conda/environment.yml
@@ -6,6 +6,7 @@ channels:
 dependencies:
   - gensim >=3.8.0
   - matplotlib
+  - numba >=0.47
   - numpy
   - pyteomics >=4.2
   - python >=3.7,<3.8
diff --git a/conda/meta.yaml b/conda/meta.yaml
index 0a4faac2..6257d9c8 100644
--- a/conda/meta.yaml
+++ b/conda/meta.yaml
@@ -36,6 +36,7 @@ requirements:
   run:
     - gensim >=3.8.0
     - matplotlib
+    - numba >=0.47
     - numpy
     - pip
     - pyteomics >=4.2
diff --git a/matchms/similarity/ModifiedCosine.py b/matchms/similarity/ModifiedCosine.py
new file mode 100644
index 00000000..e758c1b9
--- /dev/null
+++ b/matchms/similarity/ModifiedCosine.py
@@ -0,0 +1,79 @@
+from typing import Tuple
+import numpy
+from matchms.similarity.collect_peak_pairs import collect_peak_pairs
+from matchms.typing import SpectrumType
+
+
+class ModifiedCosine:
+    """Calculate 'modified cosine score' between mass spectra.
+
+    The modified cosine score aims at quantifying the similarity between two
+    mass spectra. The score is calculated by finding best possible matches between
+    peaks of two spectra. Two peaks are considered a potential match if their
+    m/z ratios lie within the given 'tolerance', or if their m/z ratios
+    lie within the tolerance once a mass-shift is applied. The mass shift is
+    simply the difference in precursor-m/z between the two spectra.
+    See Watrous et al. [PNAS, 2012, https://www.pnas.org/content/109/26/E1743]
+    for further details.
+
+    Args:
+    -----
+    tolerance: float
+        Peaks will be considered a match when <= tolerance a part. Default is 0.1.
+    """
+    def __init__(self, tolerance=0.1):
+        self.tolerance = tolerance
+
+    def __call__(self, spectrum1: SpectrumType, spectrum2: SpectrumType) -> Tuple[float, int]:
+        """Calculate modified cosine score between two spectra.
+        Args:
+        -----
+        spectrum1: SpectrumType
+            Input spectrum 1.
+        spectrum2: SpectrumType
+            Input spectrum 2.
+
+        Returns:
+        --------
+
+        Tuple with cosine score and number of matched peaks.
+        """
+        def get_peaks_arrays():
+            """Get peaks mz and intensities as numpy array."""
+            spec1 = numpy.vstack((spectrum1.peaks.mz, spectrum1.peaks.intensities)).T
+            spec2 = numpy.vstack((spectrum2.peaks.mz, spectrum2.peaks.intensities)).T
+            assert max(spec1[:, 1]) <= 1, ("Input spectrum1 is not normalized. ",
+                                           "Apply 'normalize_intensities' filter first.")
+            assert max(spec2[:, 1]) <= 1, ("Input spectrum2 is not normalized. ",
+                                           "Apply 'normalize_intensities' filter first.")
+            return spec1, spec2
+
+        def get_matching_pairs():
+            """Find all pairs of peaks that match within the given tolerance."""
+            zero_pairs = collect_peak_pairs(spec1, spec2, self.tolerance, shift=0.0)
+            message = "Precursor_mz missing. Apply 'add_precursor_mz' filter first."
+            assert spectrum1.get("precursor_mz") and spectrum2.get("precursor_mz"), message
+            mass_shift = spectrum1.get("precursor_mz") - spectrum2.get("precursor_mz")
+            nonzero_pairs = collect_peak_pairs(spec1, spec2, self.tolerance, shift=mass_shift)
+            unsorted_matching_pairs = zero_pairs + nonzero_pairs
+            return sorted(unsorted_matching_pairs, key=lambda x: x[2], reverse=True)
+
+        def calc_score():
+            """Calculate modified cosine score."""
+            used1 = set()
+            used2 = set()
+            score = 0.0
+            used_matches = []
+            for match in matching_pairs:
+                if not match[0] in used1 and not match[1] in used2:
+                    score += match[2]
+                    used1.add(match[0])  # Every peak can only be paired once
+                    used2.add(match[1])  # Every peak can only be paired once
+                    used_matches.append(match)
+            # Normalize score:
+            score = score/max(numpy.sum(spec1[:, 1]**2), numpy.sum(spec2[:, 1]**2))
+            return score, len(used_matches)
+
+        spec1, spec2 = get_peaks_arrays()
+        matching_pairs = get_matching_pairs()
+        return calc_score()
diff --git a/matchms/similarity/__init__.py b/matchms/similarity/__init__.py
index dfd559d1..a0916968 100644
--- a/matchms/similarity/__init__.py
+++ b/matchms/similarity/__init__.py
@@ -1,9 +1,11 @@
 """similarity module"""
 from .CosineGreedy import CosineGreedy
 from .IntersectMz import IntersectMz
+from .ModifiedCosine import ModifiedCosine
 
 
 __all__ = [
     "CosineGreedy",
-    "IntersectMz"
+    "IntersectMz",
+    "ModifiedCosine",
 ]
diff --git a/matchms/similarity/collect_peak_pairs.py b/matchms/similarity/collect_peak_pairs.py
new file mode 100644
index 00000000..9ba68a18
--- /dev/null
+++ b/matchms/similarity/collect_peak_pairs.py
@@ -0,0 +1,33 @@
+import numba
+import numpy
+
+
+@numba.njit
+def collect_peak_pairs(spec1, spec2, tolerance, shift=0):
+    """Find matching pairs between two spectra.
+
+    Args
+    ----
+    spec1: numpy array
+        Spectrum peaks and intensities as numpy array.
+    spec2: numpy array
+        Spectrum peaks and intensities as numpy array.
+    tolerance : float
+        Peaks will be considered a match when <= tolerance a part.
+    shift : float, optional
+        Shift spectra peaks by shift. The default is 0.
+
+    Returns
+    -------
+    matching_pairs : list
+        List of found matching peaks.
+    """
+    matching_pairs = []
+
+    for idx in range(len(spec1)):
+        intensity = spec1[idx, 1]
+        matches = numpy.where((numpy.abs(spec2[:, 0] - spec1[idx, 0] + shift) <= tolerance))[0]
+        for match in matches:
+            matching_pairs.append((idx, match, intensity*spec2[match][1]))
+
+    return matching_pairs


Test patch:
diff --git a/tests/test_collect_peak_pairs.py b/tests/test_collect_peak_pairs.py
new file mode 100644
index 00000000..d259cbe0
--- /dev/null
+++ b/tests/test_collect_peak_pairs.py
@@ -0,0 +1,28 @@
+import numpy
+from matchms.similarity.collect_peak_pairs import collect_peak_pairs
+
+
+def test_cosine_hungarian_tolerance_01():
+    """Test finding expected peak matches within tolerance=0.2."""
+    spec1 = numpy.array([[100, 200, 300, 500],
+                         [0.1, 0.1, 1.0, 1.0]], dtype="float").T
+
+    spec2 = numpy.array([[105, 205.1, 300, 500.1],
+                         [0.1, 0.1, 1.0, 1.0]], dtype="float").T
+
+    matching_pairs = collect_peak_pairs(spec1, spec2, tolerance=0.2)
+    assert len(matching_pairs) == 2, "Expected different number of matching peaks"
+    assert matching_pairs == [(2, 2, 1.0), (3, 3, 1.0)], "Expected different matchin pairs."
+
+
+def test_cosine_hungarian_tolerance_01_shift_min5():
+    """Test finding expected peak matches when given a mass_shift of -5.0."""
+    spec1 = numpy.array([[100, 200, 300, 500],
+                         [1.0, 1.0, 0.1, 0.1]], dtype="float").T
+
+    spec2 = numpy.array([[105, 205.1, 300, 500.1],
+                         [1.0, 1.0, 0.1, 0.1]], dtype="float").T
+
+    matching_pairs = collect_peak_pairs(spec1, spec2, tolerance=0.2, shift=-5.0)
+    assert len(matching_pairs) == 2, "Expected different number of matching peaks"
+    assert matching_pairs == [(0, 0, 1.0), (1, 1, 1.0)], "Expected different matchin pairs."
diff --git a/tests/test_modified_cosine.py b/tests/test_modified_cosine.py
new file mode 100644
index 00000000..bbaed356
--- /dev/null
+++ b/tests/test_modified_cosine.py
@@ -0,0 +1,82 @@
+import numpy
+import pytest
+from matchms import Spectrum
+from matchms.filtering import normalize_intensities
+from matchms.similarity import ModifiedCosine
+
+
+def test_modified_cosine_without_precursor_mz():
+    """Test without precursor-m/z. Should raise assertion error."""
+    spectrum_1 = Spectrum(mz=numpy.array([100, 150, 200, 300, 500, 510, 1100], dtype="float"),
+                          intensities=numpy.array([700, 200, 100, 1000, 200, 5, 500], dtype="float"))
+
+    spectrum_2 = Spectrum(mz=numpy.array([100, 140, 190, 300, 490, 510, 1090], dtype="float"),
+                          intensities=numpy.array([700, 200, 100, 1000, 200, 5, 500], dtype="float"))
+
+    norm_spectrum_1 = normalize_intensities(spectrum_1)
+    norm_spectrum_2 = normalize_intensities(spectrum_2)
+    modified_cosine = ModifiedCosine()
+
+    with pytest.raises(AssertionError) as msg:
+        modified_cosine(norm_spectrum_1, norm_spectrum_2)
+
+    expected_message = "Precursor_mz missing. Apply 'add_precursor_mz' filter first."
+    assert str(msg.value) == expected_message
+
+
+def test_modified_cosine_with_mass_shift_5():
+    """Test modified cosine on two spectra with mass set shift."""
+    spectrum_1 = Spectrum(mz=numpy.array([100, 150, 200, 300, 500, 510, 1100], dtype="float"),
+                          intensities=numpy.array([700, 200, 100, 1000, 200, 5, 500], dtype="float"),
+                          metadata={"precursor_mz": 1000.0})
+
+    spectrum_2 = Spectrum(mz=numpy.array([55, 105, 205, 304.5, 494.5, 515.5, 1045], dtype="float"),
+                          intensities=numpy.array([700, 200, 100, 1000, 200, 5, 500], dtype="float"),
+                          metadata={"precursor_mz": 1005.0})
+
+    norm_spectrum_1 = normalize_intensities(spectrum_1)
+    norm_spectrum_2 = normalize_intensities(spectrum_2)
+    modified_cosine = ModifiedCosine()
+    score, n_matches = modified_cosine(norm_spectrum_1, norm_spectrum_2)
+
+    assert score == pytest.approx(0.081966, 0.0001), "Expected different cosine score."
+    assert n_matches == 2, "Expected 2 matching peaks."
+
+
+def test_modified_cosine_with_mass_shift_5_tolerance_2():
+    """Test modified cosine on two spectra with mass set shift and tolerance."""
+    spectrum_1 = Spectrum(mz=numpy.array([100, 200, 299, 300, 301, 500, 510], dtype="float"),
+                          intensities=numpy.array([10, 10, 500, 100, 200, 20, 100], dtype="float"),
+                          metadata={"precursor_mz": 1000.0})
+
+    spectrum_2 = Spectrum(mz=numpy.array([105, 205, 305, 306, 505, 517], dtype="float"),
+                          intensities=numpy.array([10, 10, 500, 100, 20, 100], dtype="float"),
+                          metadata={"precursor_mz": 1005})
+
+    norm_spectrum_1 = normalize_intensities(spectrum_1)
+    norm_spectrum_2 = normalize_intensities(spectrum_2)
+    modified_cosine = ModifiedCosine(tolerance=2.0)
+    score, n_matches = modified_cosine(norm_spectrum_1, norm_spectrum_2)
+
+    assert score == pytest.approx(0.903412, 0.0001), "Expected different modified cosine score."
+    assert n_matches == 6, "Expected 6 matching peaks."
+
+
+def test_modified_cosine_order_of_input_spectrums():
+    """Test modified cosine on two spectra in changing order."""
+    spectrum_1 = Spectrum(mz=numpy.array([100, 150, 200, 300, 500, 510, 1100], dtype="float"),
+                          intensities=numpy.array([700, 200, 100, 1000, 200, 5, 500], dtype="float"),
+                          metadata={"precursor_mz": 1000.0})
+
+    spectrum_2 = Spectrum(mz=numpy.array([55, 105, 205, 304.5, 494.5, 515.5, 1045], dtype="float"),
+                          intensities=numpy.array([700, 200, 100, 1000, 200, 5, 500], dtype="float"),
+                          metadata={"precursor_mz": 1005.0})
+
+    norm_spectrum_1 = normalize_intensities(spectrum_1)
+    norm_spectrum_2 = normalize_intensities(spectrum_2)
+    modified_cosine = ModifiedCosine(tolerance=2.0)
+    score_1_2, n_matches_1_2 = modified_cosine(norm_spectrum_1, norm_spectrum_2)
+    score_2_1, n_matches_2_1 = modified_cosine(norm_spectrum_2, norm_spectrum_1)
+
+    assert score_1_2 == score_2_1, "Expected that the order of the arguments would not matter."
+    assert n_matches_1_2 == n_matches_2_1, "Expected that the order of the arguments would not matter."