corems.molecular_id.calc.SpectralSimilarity

  1__author__ = "Yuri E. Corilo"
  2__date__ = "Jun 09, 2021"
  3
  4from numpy.fft import rfft
  5from scipy.stats import pearsonr, spearmanr, kendalltau
  6from numpy import (
  7    power,
  8    dot,
  9    absolute,
 10    sqrt,
 11)
 12from numpy import sum as np_sum
 13from numpy.linalg import norm
 14from pandas import DataFrame
 15import numpy as np
 16
 17methods_name = {
 18    # "entropy_distance": "Entropy Distance",
 19    # "weighted_entropy_distance": "Dynamic weighted entropy Distance",
 20    "chebyshev_distance": "Chebyshev Distance",
 21    "squared_euclidean_distance": "Squared Euclidean Distance",
 22    "fidelity_similarity": "Fidelity Similarity",
 23    "matusita_distance": "Matusita Distance",
 24    "squared_chord_distance": "Squared-chord Distance",
 25    # "bhattacharya_1_distance": "Bhattacharya 1 Distance",
 26    # "bhattacharya_2_distance": "Bhattacharya 2 Distance",
 27    "harmonic_mean_similarity": "Harmonic mean Distance",
 28    "Pearson_chi_squared_distance": "Pearson Chi Squared Distance",
 29    "Neyman_chi_squared_distance": "Neyman Chi Squared Distance",
 30    "probabilistic_symmetric_chi_squared_distance": "Probabilistic symmetric X2 Distance",
 31    "topsoe_distance": "Topsoe Distance",
 32    "chernoff_distance": "Chernoff Distance",
 33    "ruzicka_distance": "Ruzicka Distance",
 34    "roberts_distance": "Roberts Distance",
 35    # "intersection_distance": "Intersection Distance",
 36    "motyka_distance": "Motyka Distance",
 37    "canberra_distance": "Canberra Distance",
 38    "canberra_metric": "Canberra Metric",
 39    "kulczynski_1_distance": "Kulczynski 1 Distance",
 40    # "baroni_urbani_buser_distance": "Baroni-Urbani-Buser Distance",
 41    # "penrose_size_distance": "Penrose size Distance",
 42    # "mean_character_distance": "Mean character Distance",
 43    "lorentzian_distance": "Lorentzian Distance",
 44    # "penrose_shape_distance": "Penrose shape Distance",
 45    "clark_distance": "Clark Distance",
 46    "hellinger_distance": "Hellinger Distance",
 47    "whittaker_index_of_association_distance": "Whittaker index of association Distance",
 48    # "similarity_index_distance": "Similarity Index Distance",
 49    # "improved_similarity_distance": "Improved Similarity",
 50    # "absolute_value_distance": "Absolute Value Distance",
 51    "spectral_contrast_angle_distance": "Spectral Contrast Angle",
 52    "wave_hedges_distance": "Wave Hedges Distance",
 53    "dice_similarity": "Dice Similarity",
 54    "inner_product_distance": "Inner Product Distance",
 55    "divergence_distance": "Divergence Distance",
 56    "jensen_difference_distance": "Jensen Differences Distance",
 57    "kumar_johnson_distance": "Kumar Johnson Distance",
 58    "avg_l_distance": "Avg (L1, L8) Distance",
 59    "vicis_wave_hadges_distance": "Vicis Wave Hadges Distance",
 60    "vicis_symmetric_chi_squared_1_distance": "Vicis-Symmetric X2 1 Distance",
 61    "vicis_symmetric_chi_squared_2_distance": "Vicis-Symmetric X2 2 Distance",
 62    "vicis_symmetric_chi_squared_3_distance": "Vicis-Symmetric X2 3 Distance",
 63    "max_symmetric_chi_squared_distance": "Max Symmetric Chi Squared Distance",
 64    "min_symmetric_chi_squared_distance": "Min Symmetric Chi Squared Distance",
 65    # "ms_for_id_v1": "MSforID Distance version 1",
 66    # "ms_for_id": "MSforID Distance",
 67    "additive_sym_chi_sq": "Additive Symmetric Chi Squared",
 68    "bhattacharya_distance": "Battacharya Distance",
 69    "generalized_ochiai_index": "Generalized Ochiai Index",
 70    "gower_distance": "Gower Distance",
 71    "impr_sqrt_cosine_sim": "Improved Square Root Cosine Similarity",
 72    "intersection_sim": "Intersection Similarity",
 73    "j_divergence": "J Divergence",
 74    "jensen_shannon_index": "Jensen Shannon Index",
 75    "k_divergence": "K Divergence",
 76    "VW6": "VW6",
 77    "VW5": "VW5",
 78    "VW4": "VW4",
 79    "VW3": "VW3",
 80    "VW2": "VW2",
 81    "VW1": "VW1",
 82    "taneja_divergence": "Taneja Divergence",
 83    "symmetric_chi_squared_distance": "Symmetric Chi Squared Distance",
 84    "squared_chi_squared_distance": "Squared Chi Squared Distance",
 85    "square_root_cosine_correlation": "Square Root Cosine Correlation",
 86    "sorensen_distance": "Sorensen Distance",
 87    "Minokowski_3": "Minokowski 3 Distance",
 88    "Minokowski_4": "Minokowski 4 Distance",
 89    "kumarjohnson_divergence": "Kumar Johnson Divergence",
 90    "kumarhassebrook_similarity": "Kumar Hassebrook Similarity",
 91    "kullbackleibler_divergence": "Kullback Leibler Divergence",
 92    "soergel_distance": "Soergel Distance",
 93}
 94
 95methods_scale = {
 96    "entropy": [0, np.log(4)],
 97    "weighted_entropy": [0, np.log(4)],
 98    "absolute_value": [0, 2],
 99    "avg_l": [0, 1.5],
100    "bhattacharya_1": [0, np.arccos(0) ** 2],
101    "bhattacharya_2": [0, np.inf],
102    "canberra": [0, np.inf],
103    "clark": [0, np.inf],
104    "divergence": [0, np.inf],
105    "euclidean": [0, np.sqrt(2)],
106    "hellinger": [0, np.inf],
107    "improved_similarity": [0, np.inf],
108    "lorentzian": [0, np.inf],
109    "manhattan": [0, 2],
110    "matusita": [0, np.sqrt(2)],
111    "mean_character": [0, 2],
112    "motyka": [-0.5, 0],
113    "ms_for_id": [-np.inf, 0],
114    "ms_for_id_v1": [0, np.inf],
115    "pearson_correlation": [-1, 1],
116    "penrose_shape": [0, np.sqrt(2)],
117    "penrose_size": [0, np.inf],
118    "probabilistic_symmetric_chi_squared": [0, 1],
119    "similarity_index": [0, np.inf],
120    "squared_chord": [0, 2],
121    "squared_euclidean": [0, 2],
122    "symmetric_chi_squared": [0, 0.5 * np.sqrt(2)],
123    "topsoe": [0, np.sqrt(2)],
124    "vicis_symmetric_chi_squared_3": [0, 2],
125    "wave_hedges": [0, np.inf],
126    "whittaker_index_of_association": [0, np.inf],
127}
128
129
130class SpectralSimilarity:
131    """Class containing methods for calculating spectral similarity between two mass spectra.
132
133    Parameters
134    ----------
135    ms_mz_abun_dict : dict
136        Dictionary of mass to abundance values for the experimental mass spectrum.
137    ref_obj : dict
138        Dictionary of mass to abundance values for the reference mass spectrum.
139    norm_func : function
140        Function to normalize the abundance values.
141
142    Attributes
143    ----------
144    normalize_func : function
145        Function to normalize the abundance values.
146    ms_mz_abun_dict : dict
147        Dictionary of mass to abundance values for the experimental mass spectrum.
148    ref_obj : dict
149        Dictionary of mass to abundance values for the reference mass spectrum.
150    exp_abun : list
151        List of abundance values for the experimental mass spectrum.
152    exp_mz : list
153        List of mass values for the experimental mass spectrum.
154    ref_mz : list
155        List of mass values for the reference mass spectrum.
156    ref_abun : list
157        List of abundance values for the reference mass spectrum.
158    ref_mz_abun_dict : dict
159        Dictionary of mass to abundance values for the reference mass spectrum.
160    df : DataFrame
161        DataFrame containing the experimental and reference mass spectrum data.
162    zero_filled_u_l : tuple
163        Tuple containing the experimental and reference mass spectrum data after zero filling and normalization.
164    common_mz_values : list
165        List of common mass values between the experimental and reference mass spectra.
166    n_x_y : int
167        Number of common mass values between the experimental and reference mass spectra.
168
169    Methods
170    -------
171    * nan_fill(df, fill_with=0).
172        Fill missing mass values with a given value.
173    * normalize(x, y, norm_func=sum).
174        Normalize the abundance values.
175    * weighted_cosine_correlation(a=0.5, b=1.3, nanfill=1e-10).
176        Calculate the weighted cosine correlation between the experimental and reference mass spectra.
177    * cosine_correlation().
178        Calculate the cosine correlation between the experimental and reference mass spectra.
179    * stein_scott().
180        Calculate the Stein-Scott similarity between the experimental and reference mass spectra.
181    * pearson_correlation().
182        Calculate the Pearson correlation between the experimental and reference mass spectra.
183    * spearman_correlation().
184        Calculate the Spearman correlation between the experimental and reference mass spectra.
185
186
187    """
188
189    def __init__(self, ms_mz_abun_dict, ref_obj, norm_func=sum):
190        self.normalize_func = norm_func
191        self.ms_mz_abun_dict = ms_mz_abun_dict
192        self.ref_obj = ref_obj
193
194        self.exp_abun = list(self.ms_mz_abun_dict.values())
195        self.exp_mz = list(self.ms_mz_abun_dict.keys())
196
197        self.ref_mz = self.ref_obj.get("mz")
198        self.ref_abun = self.ref_obj.get("abundance")
199
200        self.ref_mz_abun_dict = dict(zip(self.ref_mz, self.ref_abun))
201
202        # parse to dataframe, easier to zerofill and tranpose
203        self.df = DataFrame([self.ms_mz_abun_dict, self.ref_mz_abun_dict])
204
205        # fill missing mz with abundance 0
206        x, y = self.nan_fill(self.df, fill_with=1e-10)
207
208        self.zero_filled_u_l = self.normalize(x, y, norm_func=self.normalize_func)
209
210        # filter out the mass values that have zero intensities in self.exp_abun
211        exp_mz_filtered = set([k for k in self.exp_mz if self.ms_mz_abun_dict[k] != 0])
212
213        # filter out the mass values that have zero intensities in self.ref_mz
214        ref_mz_filtered = set([k for k in self.ref_mz if self.ref_mz_abun_dict[k] != 0])
215
216        # find the intersection/common mass values of both ref and exp, and sort them
217        self.common_mz_values = sorted(
218            list(exp_mz_filtered.intersection(ref_mz_filtered))
219        )
220
221        # find the number of common mass values (after filtering 0s)
222        self.n_x_y = len(self.common_mz_values)
223        # print(self.n_x_y)
224
225    def nan_fill(self, df, fill_with=0):
226        """Fill missing mass values with a given value.
227
228        Parameters
229        ----------
230        df : DataFrame
231            DataFrame containing the experimental and reference mass spectrum data.
232        fill_with : float
233            Value to fill missing mass values with.
234
235        Returns
236        -------
237        x : list
238            List of abundance values for the experimental mass spectrum.
239        y : list
240            List of abundance values for the reference mass spectrum."""
241        # CoW-safe: avoid inplace fillna; prefer to_numpy over .values.
242        # df is shape (2, n_mz) with row 0 = experimental, row 1 = reference.
243        # Equivalent to the historical `df.T[0].values, df.T[1].values`.
244        df = df.fillna(fill_with)
245        return df.iloc[0].to_numpy(), df.iloc[1].to_numpy()
246
247    def normalize(self, x, y, norm_func=sum):
248        """Normalize the abundance values.
249
250        Parameters
251        ----------
252        x : list
253            List of abundance values for the experimental mass spectrum.
254        y : list
255            List of abundance values for the reference mass spectrum.
256        norm_func : function
257            Function to normalize the abundance values.
258            Default is sum
259
260        Returns
261        -------
262        u_l : tuple
263            Tuple containing the experimental and reference mass spectrum data after zero filling and normalization.
264        """
265        if norm_func is not None:
266            u_l = (x / norm_func(x), y / norm_func(y))
267            return u_l
268        else:
269            return (x, y)
270
271    def weighted_cosine_correlation(self, a=0.5, b=1.3, nanfill=1e-10):
272        """Calculate the weighted cosine correlation between the experimental and reference mass spectra.
273
274        Parameters
275        ----------
276        a : float
277            Weighting factor for the abundance values.
278            Default is 0.5
279        b : float
280            Weighting factor for the mass values.
281            Default is 1.3
282        nanfill : float
283            Value to fill missing mass values with.
284            Default is 1e-10
285
286        Returns
287        -------
288        correlation : float
289            Weighted cosine correlation between the experimental and reference mass spectra.
290        """
291        # create dict['mz'] = abundance, for experimental data
292        # ms_mz_abun_dict = mass_spec.mz_abun_dict
293        # weight exp data
294
295        xc = power(self.exp_abun, a) * power(self.exp_mz, b)
296
297        # track back to individual mz
298        weighted_exp_dict = dict(zip(self.ms_mz_abun_dict.keys(), xc))
299
300        # weight ref data
301        yc = power(self.ref_obj.get("abundance"), a) * power(self.ref_obj.get("mz"), b)
302
303        ref_mz_abun_dict = dict(zip(self.ref_obj.get("mz"), yc))
304
305        # parse to dataframe, easier to zerofill and tranpose
306        df = DataFrame([weighted_exp_dict, ref_mz_abun_dict])
307
308        # fill missing mz with weight {abun**a}{m/z**b} to 0
309        x, y = self.nan_fill(df, fill_with=nanfill)
310
311        # correlation = (1 - cosine(x, y))
312
313        correlation = dot(x, y) / (norm(x) * norm(y))
314
315        return correlation
316
317    def cosine_correlation(self):
318        """Calculate the cosine correlation between the experimental and reference mass spectra.
319
320        Returns
321        -------
322        correlation : float
323            Cosine correlation between the experimental and reference mass spectra.
324
325        """
326        # calculate cosine correlation,
327        x = self.zero_filled_u_l[0]
328        y = self.zero_filled_u_l[1]
329
330        # correlation = (1 - cosine(x, y))
331
332        correlation = dot(x, y) / (norm(x) * norm(y))
333
334        return correlation
335
336    def stein_scott(self):
337        """Calculate the Stein-Scott similarity between the experimental and reference mass spectra.
338
339        Returns
340        -------
341        s_ss_x_y : float
342            Stein-Scott similarity between the experimental and reference mass spectra.
343        s_ss_x_y_nist : float
344            Stein-Scott similarity between the experimental and reference mass spectra.
345        """
346        # TODO check this code
347        if self.n_x_y == 0:
348            return 0, 0
349
350        # count number of non-zero abundance/peak intensity values
351        n_x = sum(a != 0 for a in self.exp_abun)
352
353        s_r_x_y = 0
354
355        a, b = 1, 0
356
357        for i in range(1, self.n_x_y):
358            current_value = self.common_mz_values[i]
359            previous_value = self.common_mz_values[i - 1]
360
361            y_i = self.ref_mz_abun_dict[current_value]
362            y_i_minus1 = self.ref_mz_abun_dict[previous_value]
363
364            lc_current = power(y_i, a) * power(current_value, b)
365            lc_previous = power(y_i_minus1, a) * power(previous_value, b)
366
367            x_i = self.ms_mz_abun_dict[current_value]
368            x_i_minus1 = self.ms_mz_abun_dict[previous_value]
369
370            uc_current = power(x_i, a) * power(current_value, b)
371            uc_previous = power(x_i_minus1, a) * power(previous_value, b)
372
373            T1 = lc_current / lc_previous
374
375            T2 = uc_previous / uc_current
376
377            temp_computation = T1 * T2
378
379            n = 0
380            if temp_computation <= 1:
381                n = 1
382            else:
383                n = -1
384
385            s_r_x_y = s_r_x_y + power(temp_computation, n)
386
387        # finish the calculation of S_R(X,Y)
388
389        s_r_x_y = s_r_x_y / self.n_x_y
390        # using the existing weighted_cosine_correlation function to get S_WC(X,Y)
391        s_wc_x_y = self.weighted_cosine_correlation(a=0.5, b=3, nanfill=0)
392
393        s_ss_x_y = ((n_x * s_wc_x_y) + (self.n_x_y * s_r_x_y)) / (n_x + self.n_x_y)
394
395        s_wc_x_y_nist = self.weighted_cosine_correlation(a=0.5, b=1.3, nanfill=0)
396
397        s_ss_x_y_nist = ((n_x * s_wc_x_y_nist) + (self.n_x_y * s_r_x_y)) / (
398            n_x + self.n_x_y
399        )
400        # final step
401
402        return s_ss_x_y, s_ss_x_y_nist
403
404    def pearson_correlation(
405        self,
406    ):
407        """Calculate the Pearson correlation between the experimental and reference mass spectra.
408
409        Returns
410        -------
411        correlation : float
412            Pearson correlation between the experimental and reference mass spectra.
413        """
414        correlation = pearsonr(self.zero_filled_u_l[0], self.zero_filled_u_l[1])
415
416        return correlation[0]
417
418    def spearman_correlation(self):
419        """Calculate the Spearman correlation between the experimental and reference mass spectra.
420
421        Returns
422        -------
423        coorelation : float
424            Spearman correlation between the experimental and reference mass spectra.
425        """
426        # calculate Spearman correlation
427        # ## TODO - Check axis
428        correlation = spearmanr(
429            self.zero_filled_u_l[0], self.zero_filled_u_l[1], axis=0
430        )
431
432        return correlation[0]
433
434    def kendall_tau(self):
435        """Calculate the Kendall's tau correlation between the experimental and reference mass spectra.
436
437        Returns
438        -------
439        correlation : float
440            Kendall's tau correlation between the experimental and reference mass spectra."""
441        # create dict['mz'] = abundance, for experimental data
442        # self.ms_mz_abun_dict = mass_spec.mz_abun_dict
443
444        # create dict['mz'] = abundance, for experimental data
445
446        # calculate Kendall's tau
447        correlation = kendalltau(self.zero_filled_u_l[0], self.zero_filled_u_l[1])
448
449        return correlation[0]
450
451    def dft_correlation(self):
452        """Calculate the DFT correlation between the experimental and reference mass spectra.
453
454        Returns
455        -------
456        correlation : float
457            DFT correlation between the experimental and reference mass spectra.
458        """
459        if self.n_x_y == 0:
460            return 0
461
462        # count number of non-zero abundance/peak intensity values
463        n_x = sum(a != 0 for a in self.exp_abun)
464
465        x, y = self.nan_fill(self.df, fill_with=0)
466
467        x, y = self.normalize(x, y, norm_func=self.normalize_func)
468
469        # get the Fourier transform of x and y
470        x_dft = rfft(x).real
471        y_dft = rfft(y).real
472
473        s_dft_xy = dot(x_dft, y_dft) / (norm(x_dft) * norm(y_dft))
474
475        # using the existing weighted_cosine_correlation function to get S_WC(X,Y)
476        s_wc_x_y = self.weighted_cosine_correlation(nanfill=0)
477
478        # final step
479        s_dft = (n_x * s_wc_x_y + self.n_x_y * s_dft_xy) / (n_x + self.n_x_y)
480
481        return s_dft
482
483    def dwt_correlation(self):
484        """Calculate the DWT correlation between the experimental and reference mass spectra.
485
486        Returns
487        -------
488        correlation : float
489            DWT correlation between the experimental and reference mass spectra.
490
491        Notes
492        -----
493        This function requires the PyWavelets library to be installed.
494            This is not a default requirement as this function is not widely used.
495        """
496
497        from pywt import dwt
498
499        if self.n_x_y == 0:
500            return 0
501
502        # count number of non-zero abundance/peak intensity values
503        n_x = sum(a != 0 for a in self.exp_abun)
504
505        # calculate cosine correlation,
506        x, y = self.nan_fill(self.df, fill_with=0)
507
508        x, y = self.normalize(x, y, norm_func=self.normalize_func)
509
510        # Make x and y into an array
511        x_a = list(x)
512        y_a = list(y)
513
514        # get the wavelet transform of x and y (Daubechies with a filter length of 4. Asymmetric. pywavelets function)
515        # Will only use the detail dwt (dwtDd
516        x_dwtD = dwt(x_a, "db2")[1]
517        y_dwtD = dwt(y_a, "db2")[1]
518
519        s_dwt_xy = dot(x_dwtD, y_dwtD) / (norm(x_dwtD) * norm(y_dwtD))
520
521        # using the existing weighted_cosine_correlation function to get S_WC(X,Y)
522        s_wc_x_y = self.weighted_cosine_correlation(nanfill=0)
523
524        # final step
525        s_dwt = (n_x * s_wc_x_y + self.n_x_y * s_dwt_xy) / (n_x + self.n_x_y)
526
527        return s_dwt
528
529    def euclidean_distance(self):
530        """Calculate the Euclidean distance between the experimental and reference mass spectra.
531
532        Returns
533        -------
534        correlation : float
535            Euclidean distance between the experimental and reference mass spectra.
536        """
537        # correlation = euclidean_distance_manual(self.zero_filled_u_l[0], self.zero_filled_u_l[1])
538        qlist = self.zero_filled_u_l[0]
539        rlist = self.zero_filled_u_l[1]
540
541        correlation = sqrt(np_sum(power(qlist - rlist, 2)))
542
543        return correlation
544
545    def manhattan_distance(self):
546        """Calculate the Manhattan distance between the experimental and reference mass spectra.
547
548        Returns
549        -------
550        correlation : float
551            Manhattan distance between the experimental and reference mass spectra.
552        """
553        qlist = self.zero_filled_u_l[0]
554        rlist = self.zero_filled_u_l[1]
555
556        return np_sum(absolute(qlist - rlist))
557
558    def jaccard_distance(self):
559        """Calculate the Jaccard distance between the experimental and reference mass spectra.
560
561        Returns
562        -------
563        correlation : float
564            Jaccard distance between the experimental and reference mass spectra.
565        """
566
567        def jaccard_similarity(list1, list2):
568            intersection = len(list(set(list1).intersection(list2)))
569            union = (len(list1) + len(list2)) - intersection
570            return float(intersection) / union
571
572        qlist = self.zero_filled_u_l[0]
573        rlist = self.zero_filled_u_l[1]
574
575        return np_sum(power(qlist - rlist, 2)) / (
576            np_sum(power(qlist, 2)) + np_sum(power(rlist, 2)) - np_sum(qlist * rlist)
577        )
578        # correlation = jaccard_similarity(self.zero_filled_u_l[0], self.zero_filled_u_l[1])
579        # @return correlation
580
581    def extra_distances(self):
582        """Function to calculate distances using additional metrics defined in math_distance.py
583
584        Currently, calculates all distances.
585
586        Returns
587        -------
588        dict_res : dict
589            Dictionary containing the distances between the experimental and reference mass spectra.
590
591        """
592        from corems.molecular_id.calc import math_distance
593
594        # qlist = self.zero_filled_u_l[2]
595        # rlist = self.zero_filled_u_l[3]
596
597        dict_res = {}
598
599        for method in methods_name:
600            # function_name = method + "_distance"
601            function_name = method
602            if hasattr(math_distance, function_name):
603                f = getattr(math_distance, function_name)
604
605                if function_name == "canberra_metric":
606                    x, y = self.nan_fill(self.df, fill_with=0)
607
608                    qlist, rlist = self.normalize(x, y, norm_func=self.normalize_func)
609                    # print("qlist:")
610                    # print(qlist)
611                    # print("rlist:")
612                    # print(rlist)
613
614                else:
615                    qlist = self.zero_filled_u_l[0]
616                    rlist = self.zero_filled_u_l[1]
617
618                dist = f(qlist, rlist)
619                # if method == "Minokowski_3":
620                #    print("qlist:")
621                #    print(qlist)
622                #    print("rlist")
623                #    print(rlist)
624                #    exit()
625                # if dist == np.nan or dis == np.inf:
626                # print(self.exp_abun)
627                # print(self.exp_mz)
628                # print(function_name)
629                # print(len(self.exp_abun))
630                # print(len(self.exp_mz))
631                # print(self.zero_filled_u_l[1])
632                dict_res[method] = dist
633
634        return dict_res