1# -*- coding: utf-8 -*-
2"""
3Created on Thu Aug 1 14:18:42 2024
4
5@author: Bart Steemans. Govers Lab.
6"""
7
8import pandas as pd
9import numpy as np
10from scipy import stats
11from skimage import filters
12from scipy.spatial.distance import pdist, squareform
13import logging
14
15bactoscoop_logger = logging.getLogger("logger")
16class SignalCorrelation:
17 def __init__(
18 self,
19 df,
20 channel1,
21 channel2,
22 feature,
23 method_name,
24 prepared_feature=None,
25 ):
26 self.df = df
27 self.channel1 = channel1
28 self.channel2 = channel2
29 self.feature = feature
30 self.method_name = method_name
31 self.prepared_feature = prepared_feature
32 self.method = self.get_method()
33 self.validate_inputs()
34
35 def get_method(self):
36 methods = {
37 "manders": self.manders_overlap_coefficient,
38 "pearson": self.pearson_correlation_coefficient,
39 "li_icq": self.li_icq,
40 "ratio": self.ratio,
41 "spearman": self.spearman_rank_correlation,
42 "kendall": self.kendall_tau,
43 "distance_corr": self.distance_correlation,
44 "covariance": self.covariance,
45 "n_cross_corr": self.normalized_cross_correlation,
46 "entropy_diff": self.entropy_difference,
47 "kurtosis_ratio": self.kurtosis_ratio,
48 "skewness_product": self.skewness_product,
49 "zero_crossings_diff": self.zero_crossings_difference,
50 "fft_peak_ratio": self.fft_peak_ratio,
51 "fft_energy_ratio": self.fft_energy_ratio,
52 "histogram_intersection": self.histogram_intersection,
53 "cosine_similarity": self.cosine_similarity,
54 }
55 if self.method_name in methods:
56 return methods[self.method_name]
57 else:
58 raise ValueError(f'Method "{self.method_name}" is not recognized')
59
60 def validate_inputs(self):
61 feature1 = f"{self.channel1}_{self.feature}"
62 feature2 = f"{self.channel2}_{self.feature}"
63 if feature1 not in self.df.columns or feature2 not in self.df.columns:
64 raise ValueError(
65 f"Features {feature1} and {feature2} must be in the DataFrame columns"
66 )
67
68 if not callable(self.method):
69 raise ValueError(f'The method "{self.method_name}" is not callable')
70
71 def calculate(self):
72 prepared_feature = self.prepared_feature
73 if prepared_feature is None:
74 prepared_feature = self.prepare_feature_pair(
75 self.df, self.channel1, self.channel2, self.feature
76 )
77
78 signals1 = prepared_feature["signals1"]
79 signals2 = prepared_feature["signals2"]
80 valid_mask = prepared_feature["valid_mask"]
81
82 results = [np.nan] * len(signals1)
83 for index, (signal1, signal2, is_valid) in enumerate(
84 zip(signals1, signals2, valid_mask)
85 ):
86 try:
87 if is_valid:
88 results[index] = self.method(signal1, signal2)
89 except ValueError as e:
90 bactoscoop_logger.debug(f"Encountered a ValueError: {e}")
91 continue
92
93 method_name = self.method.__name__
94 new_feature_name = (
95 f"{self.channel1}_{self.channel2}_{self.feature}_{method_name}"
96 )
97 self.df[new_feature_name] = results
98 return self.df
99
100 @classmethod
101 def prepare_feature_pair(cls, df, channel1, channel2, feature):
102 feature1 = f"{channel1}_{feature}"
103 feature2 = f"{channel2}_{feature}"
104 signals1 = df[feature1].to_numpy(copy=False)
105 signals2 = df[feature2].to_numpy(copy=False)
106 valid_mask = np.fromiter(
107 (
108 cls.is_valid_signal(signal1) and cls.is_valid_signal(signal2)
109 for signal1, signal2 in zip(signals1, signals2)
110 ),
111 dtype=bool,
112 count=len(signals1),
113 )
114 return {"signals1": signals1, "signals2": signals2, "valid_mask": valid_mask}
115
116 @staticmethod
117 def is_valid_signal(signal):
118 try:
119 signal_array = np.asarray(signal, dtype=float)
120 except (TypeError, ValueError):
121 return False
122 return bool(signal_array.size and np.isfinite(signal_array).all())
123
124 @staticmethod
125 def ratio(signal1, signal2):
126 if (
127 np.isscalar(signal1)
128 and np.isscalar(signal2)
129 and np.isfinite(signal1)
130 and np.isfinite(signal2)
131 and signal2 != 0
132 ):
133 return signal1 / signal2
134 return np.nan
135
136 @staticmethod
137 def manders_overlap_coefficient(signal1, signal2):
138 if isinstance(signal1, (list, np.ndarray)) and isinstance(
139 signal2, (list, np.ndarray)
140 ):
141 if len(signal1) > 1 and len(signal2) > 1:
142 signal1 = np.array(signal1)
143 signal2 = np.array(signal2)
144 return np.sum(signal1 * signal2) / np.sqrt(
145 np.sum(signal1**2) * np.sum(signal2**2)
146 )
147 return np.nan
148
149 @staticmethod
150 def pearson_correlation_coefficient(signal1, signal2):
151 if isinstance(signal1, (list, np.ndarray)) and isinstance(
152 signal2, (list, np.ndarray)
153 ):
154 if len(signal1) > 1 and len(signal2) > 1:
155 return stats.pearsonr(signal1, signal2)[0]
156 return np.nan
157
158 @staticmethod
159 def li_icq(signal1, signal2):
160 if isinstance(signal1, (list, np.ndarray)) and isinstance(
161 signal2, (list, np.ndarray)
162 ):
163 if len(signal1) > 1 and len(signal2) > 1:
164 signal1 = np.array(signal1)
165 signal2 = np.array(signal2)
166 mean1, mean2 = np.mean(signal1), np.mean(signal2)
167 product = ((signal1 - mean1) * (signal2 - mean2)) > 0
168 return (np.sum(product) / len(signal1) - 0.5) * 2
169 return np.nan
170
171 @staticmethod
172 def spearman_rank_correlation(signal1, signal2):
173 if isinstance(signal1, (list, np.ndarray)) and isinstance(
174 signal2, (list, np.ndarray)
175 ):
176 if len(signal1) > 1 and len(signal2) > 1:
177 return stats.spearmanr(signal1, signal2)[0]
178 return np.nan
179
180 @staticmethod
181 def kendall_tau(signal1, signal2):
182 if isinstance(signal1, (list, np.ndarray)) and isinstance(
183 signal2, (list, np.ndarray)
184 ):
185 if len(signal1) > 1 and len(signal2) > 1:
186 corr, _ = stats.kendalltau(signal1, signal2)
187 return corr
188 return np.nan
189
190 @staticmethod
191 def distance_correlation(signal1, signal2):
192 """Biased sample distance correlation; undefined profiles return NaN."""
193 if isinstance(signal1, (list, np.ndarray)) and isinstance(
194 signal2, (list, np.ndarray)
195 ):
196 if len(signal1) > 1 and len(signal2) > 1:
197 x = np.asarray(signal1, dtype=float)
198 y = np.asarray(signal2, dtype=float)
199 if (
200 x.ndim != 1
201 or y.ndim != 1
202 or len(x) != len(y)
203 or not np.isfinite(x).all()
204 or not np.isfinite(y).all()
205 ):
206 return np.nan
207
208 # Biased sample distance correlation uses separately centered
209 # pairwise distances for each signal.
210 distance_x = squareform(pdist(x[:, None]))
211 distance_y = squareform(pdist(y[:, None]))
212 centered_x = (
213 distance_x
214 - distance_x.mean(axis=0)
215 - distance_x.mean(axis=1)[:, None]
216 + distance_x.mean()
217 )
218 centered_y = (
219 distance_y
220 - distance_y.mean(axis=0)
221 - distance_y.mean(axis=1)[:, None]
222 + distance_y.mean()
223 )
224 variance_x = np.mean(centered_x**2)
225 variance_y = np.mean(centered_y**2)
226 if variance_x == 0 or variance_y == 0:
227 return np.nan
228 correlation_squared = np.mean(centered_x * centered_y) / np.sqrt(
229 variance_x * variance_y
230 )
231 return float(np.sqrt(np.clip(correlation_squared, 0, 1)))
232 return np.nan
233
234 @staticmethod
235 def covariance(signal1, signal2):
236 if isinstance(signal1, (list, np.ndarray)) and isinstance(
237 signal2, (list, np.ndarray)
238 ):
239 if len(signal1) > 1 and len(signal2) > 1:
240 return np.cov(signal1, signal2)[0, 1]
241 return np.nan
242
243 @staticmethod
244 def normalized_cross_correlation(signal1, signal2):
245 if isinstance(signal1, (list, np.ndarray)) and isinstance(
246 signal2, (list, np.ndarray)
247 ):
248 if len(signal1) > 1 and len(signal2) > 1:
249 return np.correlate(
250 signal1 - np.mean(signal1), signal2 - np.mean(signal2)
251 )[0] / (len(signal1) * np.std(signal1) * np.std(signal2))
252 return np.nan
253
254 @staticmethod
255 def entropy_difference(signal1, signal2):
256 if isinstance(signal1, (list, np.ndarray)) and isinstance(
257 signal2, (list, np.ndarray)
258 ):
259 if len(signal1) > 1 and len(signal2) > 1:
260 entropy1 = stats.entropy(signal1)
261 entropy2 = stats.entropy(signal2)
262 return np.abs(entropy1 - entropy2)
263 return np.nan
264
265 @staticmethod
266 def kurtosis_ratio(signal1, signal2):
267 if isinstance(signal1, (list, np.ndarray)) and isinstance(
268 signal2, (list, np.ndarray)
269 ):
270 if len(signal1) > 1 and len(signal2) > 1:
271 kurtosis1 = stats.kurtosis(signal1)
272 kurtosis2 = stats.kurtosis(signal2)
273 return np.abs(kurtosis1 / kurtosis2) if kurtosis2 != 0 else np.nan
274 return np.nan
275
276 @staticmethod
277 def skewness_product(signal1, signal2):
278 if isinstance(signal1, (list, np.ndarray)) and isinstance(
279 signal2, (list, np.ndarray)
280 ):
281 if len(signal1) > 1 and len(signal2) > 1:
282 skew1 = stats.skew(signal1)
283 skew2 = stats.skew(signal2)
284 return skew1 * skew2
285 return np.nan
286
287 @staticmethod
288 def zero_crossings_difference(signal1, signal2):
289 if isinstance(signal1, (list, np.ndarray)) and isinstance(
290 signal2, (list, np.ndarray)
291 ):
292 if len(signal1) > 1 and len(signal2) > 1:
293 zero_crossings1 = np.sum(np.diff(np.sign(signal1)) != 0)
294 zero_crossings2 = np.sum(np.diff(np.sign(signal2)) != 0)
295 return np.abs(zero_crossings1 - zero_crossings2)
296 return np.nan
297
298 @staticmethod
299 def fft_peak_ratio(signal1, signal2):
300 if isinstance(signal1, (list, np.ndarray)) and isinstance(
301 signal2, (list, np.ndarray)
302 ):
303 if len(signal1) > 1 and len(signal2) > 1:
304 fft1 = np.abs(np.fft.fft(signal1))
305 fft2 = np.abs(np.fft.fft(signal2))
306 peak1 = np.max(fft1[1:]) # Exclude DC component
307 peak2 = np.max(fft2[1:]) # Exclude DC component
308 return peak1 / peak2 if peak2 != 0 else np.nan
309 return np.nan
310
311 @staticmethod
312 def fft_energy_ratio(signal1, signal2):
313 if isinstance(signal1, (list, np.ndarray)) and isinstance(
314 signal2, (list, np.ndarray)
315 ):
316 if len(signal1) > 1 and len(signal2) > 1:
317 # Compute the FFT of the signals
318 fft1 = np.fft.fft(signal1)
319 fft2 = np.fft.fft(signal2)
320
321 # Compute the energy in each frequency bin
322 energy1 = np.sum(np.abs(fft1) ** 2)
323 energy2 = np.sum(np.abs(fft2) ** 2)
324
325 # Return the energy ratio
326 return energy1 / energy2 if energy2 != 0 else np.nan
327 return np.nan
328
329 @staticmethod
330 def histogram_intersection(signal1, signal2, bins=10):
331 """Intersect probability histograms using shared bins for both profiles."""
332 if isinstance(signal1, (list, np.ndarray)) and isinstance(
333 signal2, (list, np.ndarray)
334 ):
335 if len(signal1) > 1 and len(signal2) > 1:
336 x = np.asarray(signal1, dtype=float)
337 y = np.asarray(signal2, dtype=float)
338 if (
339 x.ndim != 1
340 or y.ndim != 1
341 or not np.isfinite(x).all()
342 or not np.isfinite(y).all()
343 ):
344 return np.nan
345
346 # Both histograms must describe the same intensity intervals.
347 # Normalize each by its own sample count so different profile
348 # lengths remain comparable and the score is symmetric.
349 edges = np.histogram_bin_edges(np.concatenate((x, y)), bins=bins)
350 hist1, _ = np.histogram(x, bins=edges)
351 hist2, _ = np.histogram(y, bins=edges)
352 count1 = hist1.sum()
353 count2 = hist2.sum()
354 if count1 == 0 or count2 == 0:
355 return np.nan
356 return float(np.minimum(hist1 / count1, hist2 / count2).sum())
357 return np.nan
358
359 @staticmethod
360 def cosine_similarity(signal1, signal2):
361 if isinstance(signal1, (list, np.ndarray)) and isinstance(
362 signal2, (list, np.ndarray)
363 ):
364 if len(signal1) > 1 and len(signal2) > 1:
365 return np.dot(signal1, signal2) / (
366 np.linalg.norm(signal1) * np.linalg.norm(signal2)
367 )
368 return np.nan