corems.chroma_peak.calc.ChromaPeakCalc
1import numpy as np 2from bisect import bisect_left 3from scipy.optimize import curve_fit 4 5try: 6 np.trapezoid 7except AttributeError: # numpy < 2.0 8 np.trapezoid = np.trapz 9 10 11__author__ = "Yuri E. Corilo" 12__date__ = "March 11, 2020" 13 14 15class GCPeakCalculation(object): 16 """ 17 Class for performing peak calculations in GC chromatography. 18 19 Methods 20 ------- 21 * `calc_area(self, tic: List[float], dx: float) -> None`: Calculate the area under the curve of the chromatogram. 22 * `linear_ri(self, right_ri: float, left_ri: float, left_rt: float, right_rt: float) -> float`: Calculate the retention index using linear interpolation. 23 * `calc_ri(self, rt_ri_pairs: List[Tuple[float, float]]) -> int`: Calculate the retention index based on the given retention time - retention index pairs. 24 """ 25 26 def calc_area(self, tic: list[float], dx: float) -> None: 27 """ 28 Calculate the area under the curve of the chromatogram. 29 30 Parameters 31 ---------- 32 tic : List[float] 33 The total ion current (TIC) values. 34 dx : float 35 The spacing between data points. 36 """ 37 yy = tic[self.start_scan : self.final_scan] 38 self._area = np.trapezoid(yy, dx=dx) 39 40 def linear_ri( 41 self, right_ri: float, left_ri: float, left_rt: float, right_rt: float 42 ) -> float: 43 """ 44 Calculate the retention index using linear interpolation. 45 46 Parameters 47 ---------- 48 right_ri : float 49 The retention index at the right reference point. 50 left_ri : float 51 The retention index at the left reference point. 52 left_rt : float 53 The retention time at the left reference point. 54 right_rt : float 55 The retention time at the right reference point. 56 57 Returns 58 ------- 59 float 60 The calculated retention index. 61 """ 62 return left_ri + ( 63 (right_ri - left_ri) 64 * (self.retention_time - left_rt) 65 / (right_rt - left_rt) 66 ) 67 68 def calc_ri(self, rt_ri_pairs: list[tuple[float, float]]) -> None: 69 """ 70 Calculate the retention index based on the given retention time - retention index pairs. 71 72 Parameters 73 ---------- 74 rt_ri_pairs : List[Tuple[float, float]] 75 The list of retention time - retention index pairs. 76 77 """ 78 current_rt = self.retention_time 79 80 rts = [rt_ri[0] for rt_ri in rt_ri_pairs] 81 index = bisect_left(rts, current_rt) 82 83 if index >= len(rt_ri_pairs): 84 index -= 1 85 86 current_ref = rt_ri_pairs[index] 87 88 if current_rt == current_ref[0]: 89 self._ri = current_ref[1] 90 91 else: 92 if index == 0: 93 index += 1 94 95 left_rt = rt_ri_pairs[index - 1][0] 96 left_ri = rt_ri_pairs[index - 1][1] 97 98 right_rt = rt_ri_pairs[index][0] 99 right_ri = rt_ri_pairs[index][1] 100 101 self._ri = self.linear_ri(right_ri, left_ri, left_rt, right_rt) 102 103 104class LCMSMassFeatureCalculation: 105 """Class for performing peak calculations in LC-MS mass spectrometry. 106 107 This class is intended to be used as a mixin class for the LCMSMassFeature class. 108 """ 109 110 def calc_dispersity_index(self): 111 """ 112 Calculate the dispersity index of the mass feature. 113 114 This function calculates the dispersity index of the mass feature and 115 stores the result in the `_dispersity_index` attribute. The dispersity index is calculated as the standard 116 deviation of the retention times that account for 50% of the cummulative intensity, starting from the most 117 intense point, as described in [1]. Note that this calculation is done within the integration bounds with 118 a pad according to the window factor, where the window factor is parameterized and encapsulated in the 119 parent LCMS object (or, if not available, defaults to 2.0 minutes before and after the apex 120 121 Returns 122 ------- 123 None, stores the result in the `_dispersity_index` attribute of the class and the `_normalized_dispersity_index` attribute, 124 which is the dispersity index normalized to the total time window used for the calculation (unitless, fraction of total window). 125 126 Raises 127 ------ 128 ValueError 129 If the EIC data are not available. 130 131 References 132 ---------- 133 1) Boiteau, Rene M., et al. "Relating Molecular Properties to the Persistence of Marine Dissolved 134 Organic Matter with Liquid Chromatography–Ultrahigh-Resolution Mass Spectrometry." 135 Environmental Science & Technology 58.7 (2024): 3267-3277. 136 """ 137 # Check if LCMSMassFeature has a parent LCMS object with a window factor 138 if hasattr(self, "mass_spectrum_obj"): 139 window_min = self.mass_spectrum_obj.parameters.lc_ms.dispersity_index_window 140 else: 141 window_min = 3.0 # minutes 142 143 # Check if the EIC data is available 144 if self.eic_list is None: 145 raise ValueError( 146 "EIC data are not available. Please add the EIC data first." 147 ) 148 149 # Define start and end of the window around the apex 150 apex_rt = self.retention_time 151 full_time = self._eic_data.time 152 full_eic = self._eic_data.eic 153 left_start = apex_rt - window_min 154 right_end = apex_rt + window_min 155 156 # Extract the EIC data within the defined window 157 time_mask = (full_time >= left_start) & (full_time <= right_end) 158 eic_subset = full_eic[time_mask] 159 time_subset = full_time[time_mask] 160 161 # Sort the EIC data and RT data by descending intensity 162 sorted_eic = eic_subset[eic_subset.argsort()[::-1]] 163 sorted_rt = time_subset[eic_subset.argsort()[::-1]] 164 165 # Calculate the dispersity index 166 cum_sum = np.cumsum(sorted_eic) / np.sum(sorted_eic) 167 rt_summ = sorted_rt[np.where(cum_sum < 0.5)] 168 if len(rt_summ) > 1: 169 d = np.std(rt_summ) 170 self._dispersity_index = d # minutes 171 self._normalized_dispersity_index = d / ( 172 time_subset[-1] - time_subset[0] 173 ) # unitless (fraction of total window used) 174 elif len(rt_summ) == 1: 175 self._dispersity_index = 0 176 self._normalized_dispersity_index = 0 177 178 def calc_fraction_height_width(self, fraction: float): 179 """ 180 Calculate the height width of the mass feature at a specfic fraction of the maximum intensity. 181 182 This function returns a tuple with the minimum and maximum half-height width based on scan resolution. 183 184 Parameters 185 ---------- 186 fraction : float 187 The fraction of the maximum intensity to calculate the height width. 188 For example, 0.5 will calculate the half-height width. 189 190 Returns 191 ------- 192 Tuple[float, float, bool] 193 The minimum and maximum half-height width based on scan resolution (in minutes), and a boolean indicating if the width was estimated. 194 """ 195 196 # Pull out the EIC data 197 eic = self._eic_data.eic_smoothed 198 199 # Find the indices of the maximum intensity on either side 200 max_index = np.where(self._eic_data.scans == self.apex_scan)[0][0] 201 left_index = max_index 202 right_index = max_index 203 while eic[left_index] > eic[max_index] * fraction and left_index > 0: 204 left_index -= 1 205 while ( 206 eic[right_index] > eic[max_index] * fraction and right_index < len(eic) - 1 207 ): 208 right_index += 1 209 210 # Get the retention times of the indexes just below the half height 211 left_rt = self._eic_data.time[left_index] 212 right_rt = self._eic_data.time[right_index] 213 214 # If left_rt and right_rt are outside the bounds of the integration, set them to the bounds and set estimated to True 215 estimated = False 216 if left_rt < self.eic_rt_list[0]: 217 left_rt = self.eic_rt_list[0] 218 left_index = np.where(self._eic_data.scans == self._eic_data.apexes[0][0])[ 219 0 220 ][0] 221 estimated = True 222 if right_rt > self.eic_rt_list[-1]: 223 right_rt = self.eic_rt_list[-1] 224 right_index = np.where( 225 self._eic_data.scans == self._eic_data.apexes[0][-1] 226 )[0][0] 227 estimated = True 228 half_height_width_max = right_rt - left_rt 229 230 # Get the retention times of the indexes just above the half height 231 left_rt = self._eic_data.time[left_index + 1] 232 right_rt = self._eic_data.time[right_index - 1] 233 half_height_width_min = right_rt - left_rt 234 235 return half_height_width_min, half_height_width_max, estimated 236 237 def calc_half_height_width(self, accept_estimated: bool = False): 238 """ 239 Calculate the half-height width of the mass feature. 240 241 This function calculates the half-height width of the mass feature and 242 stores the result in the `_half_height_width` attribute 243 244 Returns 245 ------- 246 None, stores the result in the `_half_height_width` attribute of the class. 247 """ 248 min_, max_, estimated = self.calc_fraction_height_width(0.5) 249 if not estimated or accept_estimated: 250 self._half_height_width = np.array([min_, max_]) 251 252 def calc_tailing_factor(self, accept_estimated: bool = False): 253 """ 254 Calculate the peak asymmetry of the mass feature. 255 256 This function calculates the peak asymmetry of the mass feature and 257 stores the result in the `_tailing_factor` attribute. 258 Calculations completed at 5% of the peak height in accordance with the USP tailing factor calculation. 259 260 Returns 261 ------- 262 None, stores the result in the `_tailing_factor` attribute of the class. 263 264 References 265 ---------- 266 1) JIS K0124:2011 General rules for high performance liquid chromatography 267 2) JIS K0214:2013 Technical terms for analytical chemistry 268 """ 269 # First calculate the width of the peak at 5% of the peak height 270 width_min, width_max, estimated = self.calc_fraction_height_width(0.05) 271 272 if not estimated or accept_estimated: 273 # Next calculate the width of the peak at 95% of the peak height 274 eic = self._eic_data.eic_smoothed 275 max_index = np.where(self._eic_data.scans == self.apex_scan)[0][0] 276 left_index = max_index 277 while eic[left_index] > eic[max_index] * 0.05 and left_index > 0: 278 left_index -= 1 279 280 left_half_time_min = ( 281 self._eic_data.time[max_index] - self._eic_data.time[left_index] 282 ) 283 left_half_time_max = ( 284 self._eic_data.time[max_index] - self._eic_data.time[left_index + 1] 285 ) 286 287 tailing_factor = np.mean([width_min, width_max]) / ( 288 2 * np.mean([left_half_time_min, left_half_time_max]) 289 ) 290 291 self._tailing_factor = tailing_factor 292 293 def calc_gaussian_similarity(self): 294 """ 295 Calculate the Gaussian similarity score of the mass feature. 296 297 This function fits a Gaussian curve to the EIC data and evaluates 298 the goodness of fit using R-squared. Note that this only uses data within 299 the set integration bounds of the mass feature. A score close to 1 indicates 300 the peak closely resembles an ideal Gaussian shape. 301 302 Returns 303 ------- 304 None, stores the result in the `_gaussian_similarity` attribute of the class. 305 306 Raises 307 ------ 308 ValueError 309 If the EIC data are not available. 310 """ 311 # Check if the EIC data is available 312 if self.eic_list is None: 313 raise ValueError( 314 "EIC data are not available. Please add the EIC data first." 315 ) 316 317 # Get EIC data within integration bounds 318 time_data = np.array(self.eic_rt_list) 319 intensity_data = np.array(self.eic_list) 320 321 if len(time_data) < 4: # Need minimum points for meaningful fit 322 self._gaussian_similarity = np.nan 323 return 324 325 # Check for valid intensity data 326 max_intensity = np.max(intensity_data) 327 if max_intensity == 0: 328 self._gaussian_similarity = np.nan 329 return 330 331 try: 332 # Define Gaussian function 333 def gaussian(x, amplitude, mean, stddev, baseline): 334 return ( 335 amplitude * np.exp(-((x - mean) ** 2) / (2 * stddev**2)) + baseline 336 ) 337 338 # Initial parameter estimates 339 amplitude_init = max_intensity 340 mean_init = time_data[np.argmax(intensity_data)] 341 stddev_init = (time_data[-1] - time_data[0]) / 6 # Rough estimate 342 baseline_init = np.min(intensity_data) 343 344 # Fit Gaussian curve 345 popt, _ = curve_fit( 346 gaussian, 347 time_data, 348 intensity_data, 349 p0=[amplitude_init, mean_init, stddev_init, baseline_init], 350 maxfev=1000, 351 bounds=( 352 [0, time_data[0], 0, 0], # Lower bounds 353 [np.inf, time_data[-1], np.inf, max_intensity], # Upper bounds 354 ), 355 ) 356 357 # Calculate fitted values 358 fitted_intensities = gaussian(time_data, *popt) 359 360 # Calculate R-squared (coefficient of determination) 361 ss_res = np.sum((intensity_data - fitted_intensities) ** 2) 362 ss_tot = np.sum((intensity_data - np.mean(intensity_data)) ** 2) 363 364 if ss_tot == 0: 365 self._gaussian_similarity = np.nan 366 else: 367 r_squared = 1 - (ss_res / ss_tot) 368 # R² should be between 0 and 1 for reasonable fits 369 # If negative, the model is worse than the mean - treat as non-computable 370 self._gaussian_similarity = r_squared if r_squared >= 0 else np.nan 371 372 except (RuntimeError, ValueError, TypeError): 373 # Fitting failed, assign NaN 374 self._gaussian_similarity = np.nan 375 376 def calc_noise_score(self): 377 """ 378 Calculate the noise score of the mass feature separately for left and right sides. 379 380 This function estimates the signal-to-noise ratio by comparing the peak 381 intensity to the baseline noise level in surrounding regions. It calculates 382 separate scores for the left and right sides of the peak, which are stored as a tuple 383 in the `_noise_score` attribute. The noise estimation windows are encapsulated in the 384 parent LCMS object (or, if not available, defaults to twice the peak width on each side). 385 386 387 Returns 388 ------- 389 None, stores the result in the `_noise_score` attribute as a tuple (left_score, right_score). 390 391 Raises 392 ------ 393 ValueError 394 If the EIC data are not available. 395 """ 396 # Check if the EIC data is available 397 if self.eic_list is None: 398 raise ValueError( 399 "EIC data are not available. Please add the EIC data first." 400 ) 401 402 # Check if LCMSMassFeature has a parent LCMS object with a window factor 403 if hasattr(self, "mass_spectrum_obj"): 404 noise_window_factor = ( 405 self.mass_spectrum_obj.parameters.lc_ms.noise_window_factor 406 ) 407 else: 408 noise_window_factor = 2.0 # times the peak width 409 410 # Get full EIC data (not just integration bounds) 411 full_time = self._eic_data.time 412 full_eic = self._eic_data.eic 413 414 # Get peak information 415 apex_rt = self.retention_time 416 peak_intensity = np.max(self.eic_list) 417 418 # Retrieve width from integration bounds 419 peak_width = self.eic_rt_list[-1] - self.eic_rt_list[0] 420 421 # Define noise estimation windows 422 noise_window_size = peak_width * noise_window_factor # in minutes 423 left_noise_start = apex_rt - peak_width - noise_window_size 424 left_noise_end = apex_rt - peak_width 425 right_noise_start = apex_rt + peak_width 426 right_noise_end = apex_rt + peak_width + noise_window_size 427 428 # Extract noise regions 429 left_noise_mask = (full_time >= left_noise_start) & ( 430 full_time <= left_noise_end 431 ) 432 right_noise_mask = (full_time >= right_noise_start) & ( 433 full_time <= right_noise_end 434 ) 435 436 left_noise = full_eic[left_noise_mask] 437 right_noise = full_eic[right_noise_mask] 438 439 # Calculate left noise score 440 if len(left_noise) == 0: 441 left_score = np.nan 442 else: 443 left_baseline = np.median(left_noise) 444 left_noise_std = np.std(left_noise) 445 446 if left_noise_std == 0: 447 if peak_intensity > left_baseline: 448 left_score = 1.0 449 else: 450 left_score = np.nan 451 else: 452 left_signal = peak_intensity - left_baseline 453 if left_signal <= 0: 454 left_score = 0.0 455 else: 456 left_snr = left_signal / left_noise_std 457 left_score = min(1.0, left_snr / (left_snr + 10.0)) 458 459 # Calculate right noise score 460 if len(right_noise) == 0: 461 right_score = np.nan 462 else: 463 right_baseline = np.median(right_noise) 464 right_noise_std = np.std(right_noise) 465 466 if right_noise_std == 0: 467 if peak_intensity > right_baseline: 468 right_score = 1.0 469 else: 470 right_score = np.nan 471 else: 472 right_signal = peak_intensity - right_baseline 473 if right_signal <= 0: 474 right_score = 0.0 475 else: 476 right_snr = right_signal / right_noise_std 477 right_score = min(1.0, right_snr / (right_snr + 10.0)) 478 479 # Store as tuple 480 self._noise_score = (left_score, right_score)
class
GCPeakCalculation:
16class GCPeakCalculation(object): 17 """ 18 Class for performing peak calculations in GC chromatography. 19 20 Methods 21 ------- 22 * `calc_area(self, tic: List[float], dx: float) -> None`: Calculate the area under the curve of the chromatogram. 23 * `linear_ri(self, right_ri: float, left_ri: float, left_rt: float, right_rt: float) -> float`: Calculate the retention index using linear interpolation. 24 * `calc_ri(self, rt_ri_pairs: List[Tuple[float, float]]) -> int`: Calculate the retention index based on the given retention time - retention index pairs. 25 """ 26 27 def calc_area(self, tic: list[float], dx: float) -> None: 28 """ 29 Calculate the area under the curve of the chromatogram. 30 31 Parameters 32 ---------- 33 tic : List[float] 34 The total ion current (TIC) values. 35 dx : float 36 The spacing between data points. 37 """ 38 yy = tic[self.start_scan : self.final_scan] 39 self._area = np.trapezoid(yy, dx=dx) 40 41 def linear_ri( 42 self, right_ri: float, left_ri: float, left_rt: float, right_rt: float 43 ) -> float: 44 """ 45 Calculate the retention index using linear interpolation. 46 47 Parameters 48 ---------- 49 right_ri : float 50 The retention index at the right reference point. 51 left_ri : float 52 The retention index at the left reference point. 53 left_rt : float 54 The retention time at the left reference point. 55 right_rt : float 56 The retention time at the right reference point. 57 58 Returns 59 ------- 60 float 61 The calculated retention index. 62 """ 63 return left_ri + ( 64 (right_ri - left_ri) 65 * (self.retention_time - left_rt) 66 / (right_rt - left_rt) 67 ) 68 69 def calc_ri(self, rt_ri_pairs: list[tuple[float, float]]) -> None: 70 """ 71 Calculate the retention index based on the given retention time - retention index pairs. 72 73 Parameters 74 ---------- 75 rt_ri_pairs : List[Tuple[float, float]] 76 The list of retention time - retention index pairs. 77 78 """ 79 current_rt = self.retention_time 80 81 rts = [rt_ri[0] for rt_ri in rt_ri_pairs] 82 index = bisect_left(rts, current_rt) 83 84 if index >= len(rt_ri_pairs): 85 index -= 1 86 87 current_ref = rt_ri_pairs[index] 88 89 if current_rt == current_ref[0]: 90 self._ri = current_ref[1] 91 92 else: 93 if index == 0: 94 index += 1 95 96 left_rt = rt_ri_pairs[index - 1][0] 97 left_ri = rt_ri_pairs[index - 1][1] 98 99 right_rt = rt_ri_pairs[index][0] 100 right_ri = rt_ri_pairs[index][1] 101 102 self._ri = self.linear_ri(right_ri, left_ri, left_rt, right_rt)
Class for performing peak calculations in GC chromatography.
Methods
calc_area(self, tic: List[float], dx: float) -> None: Calculate the area under the curve of the chromatogram.linear_ri(self, right_ri: float, left_ri: float, left_rt: float, right_rt: float) -> float: Calculate the retention index using linear interpolation.calc_ri(self, rt_ri_pairs: List[Tuple[float, float]]) -> int: Calculate the retention index based on the given retention time - retention index pairs.
def
calc_area(self, tic: list[float], dx: float) -> None:
27 def calc_area(self, tic: list[float], dx: float) -> None: 28 """ 29 Calculate the area under the curve of the chromatogram. 30 31 Parameters 32 ---------- 33 tic : List[float] 34 The total ion current (TIC) values. 35 dx : float 36 The spacing between data points. 37 """ 38 yy = tic[self.start_scan : self.final_scan] 39 self._area = np.trapezoid(yy, dx=dx)
Calculate the area under the curve of the chromatogram.
Parameters
- tic (List[float]): The total ion current (TIC) values.
- dx (float): The spacing between data points.
def
linear_ri( self, right_ri: float, left_ri: float, left_rt: float, right_rt: float) -> float:
41 def linear_ri( 42 self, right_ri: float, left_ri: float, left_rt: float, right_rt: float 43 ) -> float: 44 """ 45 Calculate the retention index using linear interpolation. 46 47 Parameters 48 ---------- 49 right_ri : float 50 The retention index at the right reference point. 51 left_ri : float 52 The retention index at the left reference point. 53 left_rt : float 54 The retention time at the left reference point. 55 right_rt : float 56 The retention time at the right reference point. 57 58 Returns 59 ------- 60 float 61 The calculated retention index. 62 """ 63 return left_ri + ( 64 (right_ri - left_ri) 65 * (self.retention_time - left_rt) 66 / (right_rt - left_rt) 67 )
Calculate the retention index using linear interpolation.
Parameters
- right_ri (float): The retention index at the right reference point.
- left_ri (float): The retention index at the left reference point.
- left_rt (float): The retention time at the left reference point.
- right_rt (float): The retention time at the right reference point.
Returns
- float: The calculated retention index.
def
calc_ri(self, rt_ri_pairs: list[tuple[float, float]]) -> None:
69 def calc_ri(self, rt_ri_pairs: list[tuple[float, float]]) -> None: 70 """ 71 Calculate the retention index based on the given retention time - retention index pairs. 72 73 Parameters 74 ---------- 75 rt_ri_pairs : List[Tuple[float, float]] 76 The list of retention time - retention index pairs. 77 78 """ 79 current_rt = self.retention_time 80 81 rts = [rt_ri[0] for rt_ri in rt_ri_pairs] 82 index = bisect_left(rts, current_rt) 83 84 if index >= len(rt_ri_pairs): 85 index -= 1 86 87 current_ref = rt_ri_pairs[index] 88 89 if current_rt == current_ref[0]: 90 self._ri = current_ref[1] 91 92 else: 93 if index == 0: 94 index += 1 95 96 left_rt = rt_ri_pairs[index - 1][0] 97 left_ri = rt_ri_pairs[index - 1][1] 98 99 right_rt = rt_ri_pairs[index][0] 100 right_ri = rt_ri_pairs[index][1] 101 102 self._ri = self.linear_ri(right_ri, left_ri, left_rt, right_rt)
Calculate the retention index based on the given retention time - retention index pairs.
Parameters
- rt_ri_pairs (List[Tuple[float, float]]): The list of retention time - retention index pairs.
class
LCMSMassFeatureCalculation:
105class LCMSMassFeatureCalculation: 106 """Class for performing peak calculations in LC-MS mass spectrometry. 107 108 This class is intended to be used as a mixin class for the LCMSMassFeature class. 109 """ 110 111 def calc_dispersity_index(self): 112 """ 113 Calculate the dispersity index of the mass feature. 114 115 This function calculates the dispersity index of the mass feature and 116 stores the result in the `_dispersity_index` attribute. The dispersity index is calculated as the standard 117 deviation of the retention times that account for 50% of the cummulative intensity, starting from the most 118 intense point, as described in [1]. Note that this calculation is done within the integration bounds with 119 a pad according to the window factor, where the window factor is parameterized and encapsulated in the 120 parent LCMS object (or, if not available, defaults to 2.0 minutes before and after the apex 121 122 Returns 123 ------- 124 None, stores the result in the `_dispersity_index` attribute of the class and the `_normalized_dispersity_index` attribute, 125 which is the dispersity index normalized to the total time window used for the calculation (unitless, fraction of total window). 126 127 Raises 128 ------ 129 ValueError 130 If the EIC data are not available. 131 132 References 133 ---------- 134 1) Boiteau, Rene M., et al. "Relating Molecular Properties to the Persistence of Marine Dissolved 135 Organic Matter with Liquid Chromatography–Ultrahigh-Resolution Mass Spectrometry." 136 Environmental Science & Technology 58.7 (2024): 3267-3277. 137 """ 138 # Check if LCMSMassFeature has a parent LCMS object with a window factor 139 if hasattr(self, "mass_spectrum_obj"): 140 window_min = self.mass_spectrum_obj.parameters.lc_ms.dispersity_index_window 141 else: 142 window_min = 3.0 # minutes 143 144 # Check if the EIC data is available 145 if self.eic_list is None: 146 raise ValueError( 147 "EIC data are not available. Please add the EIC data first." 148 ) 149 150 # Define start and end of the window around the apex 151 apex_rt = self.retention_time 152 full_time = self._eic_data.time 153 full_eic = self._eic_data.eic 154 left_start = apex_rt - window_min 155 right_end = apex_rt + window_min 156 157 # Extract the EIC data within the defined window 158 time_mask = (full_time >= left_start) & (full_time <= right_end) 159 eic_subset = full_eic[time_mask] 160 time_subset = full_time[time_mask] 161 162 # Sort the EIC data and RT data by descending intensity 163 sorted_eic = eic_subset[eic_subset.argsort()[::-1]] 164 sorted_rt = time_subset[eic_subset.argsort()[::-1]] 165 166 # Calculate the dispersity index 167 cum_sum = np.cumsum(sorted_eic) / np.sum(sorted_eic) 168 rt_summ = sorted_rt[np.where(cum_sum < 0.5)] 169 if len(rt_summ) > 1: 170 d = np.std(rt_summ) 171 self._dispersity_index = d # minutes 172 self._normalized_dispersity_index = d / ( 173 time_subset[-1] - time_subset[0] 174 ) # unitless (fraction of total window used) 175 elif len(rt_summ) == 1: 176 self._dispersity_index = 0 177 self._normalized_dispersity_index = 0 178 179 def calc_fraction_height_width(self, fraction: float): 180 """ 181 Calculate the height width of the mass feature at a specfic fraction of the maximum intensity. 182 183 This function returns a tuple with the minimum and maximum half-height width based on scan resolution. 184 185 Parameters 186 ---------- 187 fraction : float 188 The fraction of the maximum intensity to calculate the height width. 189 For example, 0.5 will calculate the half-height width. 190 191 Returns 192 ------- 193 Tuple[float, float, bool] 194 The minimum and maximum half-height width based on scan resolution (in minutes), and a boolean indicating if the width was estimated. 195 """ 196 197 # Pull out the EIC data 198 eic = self._eic_data.eic_smoothed 199 200 # Find the indices of the maximum intensity on either side 201 max_index = np.where(self._eic_data.scans == self.apex_scan)[0][0] 202 left_index = max_index 203 right_index = max_index 204 while eic[left_index] > eic[max_index] * fraction and left_index > 0: 205 left_index -= 1 206 while ( 207 eic[right_index] > eic[max_index] * fraction and right_index < len(eic) - 1 208 ): 209 right_index += 1 210 211 # Get the retention times of the indexes just below the half height 212 left_rt = self._eic_data.time[left_index] 213 right_rt = self._eic_data.time[right_index] 214 215 # If left_rt and right_rt are outside the bounds of the integration, set them to the bounds and set estimated to True 216 estimated = False 217 if left_rt < self.eic_rt_list[0]: 218 left_rt = self.eic_rt_list[0] 219 left_index = np.where(self._eic_data.scans == self._eic_data.apexes[0][0])[ 220 0 221 ][0] 222 estimated = True 223 if right_rt > self.eic_rt_list[-1]: 224 right_rt = self.eic_rt_list[-1] 225 right_index = np.where( 226 self._eic_data.scans == self._eic_data.apexes[0][-1] 227 )[0][0] 228 estimated = True 229 half_height_width_max = right_rt - left_rt 230 231 # Get the retention times of the indexes just above the half height 232 left_rt = self._eic_data.time[left_index + 1] 233 right_rt = self._eic_data.time[right_index - 1] 234 half_height_width_min = right_rt - left_rt 235 236 return half_height_width_min, half_height_width_max, estimated 237 238 def calc_half_height_width(self, accept_estimated: bool = False): 239 """ 240 Calculate the half-height width of the mass feature. 241 242 This function calculates the half-height width of the mass feature and 243 stores the result in the `_half_height_width` attribute 244 245 Returns 246 ------- 247 None, stores the result in the `_half_height_width` attribute of the class. 248 """ 249 min_, max_, estimated = self.calc_fraction_height_width(0.5) 250 if not estimated or accept_estimated: 251 self._half_height_width = np.array([min_, max_]) 252 253 def calc_tailing_factor(self, accept_estimated: bool = False): 254 """ 255 Calculate the peak asymmetry of the mass feature. 256 257 This function calculates the peak asymmetry of the mass feature and 258 stores the result in the `_tailing_factor` attribute. 259 Calculations completed at 5% of the peak height in accordance with the USP tailing factor calculation. 260 261 Returns 262 ------- 263 None, stores the result in the `_tailing_factor` attribute of the class. 264 265 References 266 ---------- 267 1) JIS K0124:2011 General rules for high performance liquid chromatography 268 2) JIS K0214:2013 Technical terms for analytical chemistry 269 """ 270 # First calculate the width of the peak at 5% of the peak height 271 width_min, width_max, estimated = self.calc_fraction_height_width(0.05) 272 273 if not estimated or accept_estimated: 274 # Next calculate the width of the peak at 95% of the peak height 275 eic = self._eic_data.eic_smoothed 276 max_index = np.where(self._eic_data.scans == self.apex_scan)[0][0] 277 left_index = max_index 278 while eic[left_index] > eic[max_index] * 0.05 and left_index > 0: 279 left_index -= 1 280 281 left_half_time_min = ( 282 self._eic_data.time[max_index] - self._eic_data.time[left_index] 283 ) 284 left_half_time_max = ( 285 self._eic_data.time[max_index] - self._eic_data.time[left_index + 1] 286 ) 287 288 tailing_factor = np.mean([width_min, width_max]) / ( 289 2 * np.mean([left_half_time_min, left_half_time_max]) 290 ) 291 292 self._tailing_factor = tailing_factor 293 294 def calc_gaussian_similarity(self): 295 """ 296 Calculate the Gaussian similarity score of the mass feature. 297 298 This function fits a Gaussian curve to the EIC data and evaluates 299 the goodness of fit using R-squared. Note that this only uses data within 300 the set integration bounds of the mass feature. A score close to 1 indicates 301 the peak closely resembles an ideal Gaussian shape. 302 303 Returns 304 ------- 305 None, stores the result in the `_gaussian_similarity` attribute of the class. 306 307 Raises 308 ------ 309 ValueError 310 If the EIC data are not available. 311 """ 312 # Check if the EIC data is available 313 if self.eic_list is None: 314 raise ValueError( 315 "EIC data are not available. Please add the EIC data first." 316 ) 317 318 # Get EIC data within integration bounds 319 time_data = np.array(self.eic_rt_list) 320 intensity_data = np.array(self.eic_list) 321 322 if len(time_data) < 4: # Need minimum points for meaningful fit 323 self._gaussian_similarity = np.nan 324 return 325 326 # Check for valid intensity data 327 max_intensity = np.max(intensity_data) 328 if max_intensity == 0: 329 self._gaussian_similarity = np.nan 330 return 331 332 try: 333 # Define Gaussian function 334 def gaussian(x, amplitude, mean, stddev, baseline): 335 return ( 336 amplitude * np.exp(-((x - mean) ** 2) / (2 * stddev**2)) + baseline 337 ) 338 339 # Initial parameter estimates 340 amplitude_init = max_intensity 341 mean_init = time_data[np.argmax(intensity_data)] 342 stddev_init = (time_data[-1] - time_data[0]) / 6 # Rough estimate 343 baseline_init = np.min(intensity_data) 344 345 # Fit Gaussian curve 346 popt, _ = curve_fit( 347 gaussian, 348 time_data, 349 intensity_data, 350 p0=[amplitude_init, mean_init, stddev_init, baseline_init], 351 maxfev=1000, 352 bounds=( 353 [0, time_data[0], 0, 0], # Lower bounds 354 [np.inf, time_data[-1], np.inf, max_intensity], # Upper bounds 355 ), 356 ) 357 358 # Calculate fitted values 359 fitted_intensities = gaussian(time_data, *popt) 360 361 # Calculate R-squared (coefficient of determination) 362 ss_res = np.sum((intensity_data - fitted_intensities) ** 2) 363 ss_tot = np.sum((intensity_data - np.mean(intensity_data)) ** 2) 364 365 if ss_tot == 0: 366 self._gaussian_similarity = np.nan 367 else: 368 r_squared = 1 - (ss_res / ss_tot) 369 # R² should be between 0 and 1 for reasonable fits 370 # If negative, the model is worse than the mean - treat as non-computable 371 self._gaussian_similarity = r_squared if r_squared >= 0 else np.nan 372 373 except (RuntimeError, ValueError, TypeError): 374 # Fitting failed, assign NaN 375 self._gaussian_similarity = np.nan 376 377 def calc_noise_score(self): 378 """ 379 Calculate the noise score of the mass feature separately for left and right sides. 380 381 This function estimates the signal-to-noise ratio by comparing the peak 382 intensity to the baseline noise level in surrounding regions. It calculates 383 separate scores for the left and right sides of the peak, which are stored as a tuple 384 in the `_noise_score` attribute. The noise estimation windows are encapsulated in the 385 parent LCMS object (or, if not available, defaults to twice the peak width on each side). 386 387 388 Returns 389 ------- 390 None, stores the result in the `_noise_score` attribute as a tuple (left_score, right_score). 391 392 Raises 393 ------ 394 ValueError 395 If the EIC data are not available. 396 """ 397 # Check if the EIC data is available 398 if self.eic_list is None: 399 raise ValueError( 400 "EIC data are not available. Please add the EIC data first." 401 ) 402 403 # Check if LCMSMassFeature has a parent LCMS object with a window factor 404 if hasattr(self, "mass_spectrum_obj"): 405 noise_window_factor = ( 406 self.mass_spectrum_obj.parameters.lc_ms.noise_window_factor 407 ) 408 else: 409 noise_window_factor = 2.0 # times the peak width 410 411 # Get full EIC data (not just integration bounds) 412 full_time = self._eic_data.time 413 full_eic = self._eic_data.eic 414 415 # Get peak information 416 apex_rt = self.retention_time 417 peak_intensity = np.max(self.eic_list) 418 419 # Retrieve width from integration bounds 420 peak_width = self.eic_rt_list[-1] - self.eic_rt_list[0] 421 422 # Define noise estimation windows 423 noise_window_size = peak_width * noise_window_factor # in minutes 424 left_noise_start = apex_rt - peak_width - noise_window_size 425 left_noise_end = apex_rt - peak_width 426 right_noise_start = apex_rt + peak_width 427 right_noise_end = apex_rt + peak_width + noise_window_size 428 429 # Extract noise regions 430 left_noise_mask = (full_time >= left_noise_start) & ( 431 full_time <= left_noise_end 432 ) 433 right_noise_mask = (full_time >= right_noise_start) & ( 434 full_time <= right_noise_end 435 ) 436 437 left_noise = full_eic[left_noise_mask] 438 right_noise = full_eic[right_noise_mask] 439 440 # Calculate left noise score 441 if len(left_noise) == 0: 442 left_score = np.nan 443 else: 444 left_baseline = np.median(left_noise) 445 left_noise_std = np.std(left_noise) 446 447 if left_noise_std == 0: 448 if peak_intensity > left_baseline: 449 left_score = 1.0 450 else: 451 left_score = np.nan 452 else: 453 left_signal = peak_intensity - left_baseline 454 if left_signal <= 0: 455 left_score = 0.0 456 else: 457 left_snr = left_signal / left_noise_std 458 left_score = min(1.0, left_snr / (left_snr + 10.0)) 459 460 # Calculate right noise score 461 if len(right_noise) == 0: 462 right_score = np.nan 463 else: 464 right_baseline = np.median(right_noise) 465 right_noise_std = np.std(right_noise) 466 467 if right_noise_std == 0: 468 if peak_intensity > right_baseline: 469 right_score = 1.0 470 else: 471 right_score = np.nan 472 else: 473 right_signal = peak_intensity - right_baseline 474 if right_signal <= 0: 475 right_score = 0.0 476 else: 477 right_snr = right_signal / right_noise_std 478 right_score = min(1.0, right_snr / (right_snr + 10.0)) 479 480 # Store as tuple 481 self._noise_score = (left_score, right_score)
Class for performing peak calculations in LC-MS mass spectrometry.
This class is intended to be used as a mixin class for the LCMSMassFeature class.
def
calc_dispersity_index(self):
111 def calc_dispersity_index(self): 112 """ 113 Calculate the dispersity index of the mass feature. 114 115 This function calculates the dispersity index of the mass feature and 116 stores the result in the `_dispersity_index` attribute. The dispersity index is calculated as the standard 117 deviation of the retention times that account for 50% of the cummulative intensity, starting from the most 118 intense point, as described in [1]. Note that this calculation is done within the integration bounds with 119 a pad according to the window factor, where the window factor is parameterized and encapsulated in the 120 parent LCMS object (or, if not available, defaults to 2.0 minutes before and after the apex 121 122 Returns 123 ------- 124 None, stores the result in the `_dispersity_index` attribute of the class and the `_normalized_dispersity_index` attribute, 125 which is the dispersity index normalized to the total time window used for the calculation (unitless, fraction of total window). 126 127 Raises 128 ------ 129 ValueError 130 If the EIC data are not available. 131 132 References 133 ---------- 134 1) Boiteau, Rene M., et al. "Relating Molecular Properties to the Persistence of Marine Dissolved 135 Organic Matter with Liquid Chromatography–Ultrahigh-Resolution Mass Spectrometry." 136 Environmental Science & Technology 58.7 (2024): 3267-3277. 137 """ 138 # Check if LCMSMassFeature has a parent LCMS object with a window factor 139 if hasattr(self, "mass_spectrum_obj"): 140 window_min = self.mass_spectrum_obj.parameters.lc_ms.dispersity_index_window 141 else: 142 window_min = 3.0 # minutes 143 144 # Check if the EIC data is available 145 if self.eic_list is None: 146 raise ValueError( 147 "EIC data are not available. Please add the EIC data first." 148 ) 149 150 # Define start and end of the window around the apex 151 apex_rt = self.retention_time 152 full_time = self._eic_data.time 153 full_eic = self._eic_data.eic 154 left_start = apex_rt - window_min 155 right_end = apex_rt + window_min 156 157 # Extract the EIC data within the defined window 158 time_mask = (full_time >= left_start) & (full_time <= right_end) 159 eic_subset = full_eic[time_mask] 160 time_subset = full_time[time_mask] 161 162 # Sort the EIC data and RT data by descending intensity 163 sorted_eic = eic_subset[eic_subset.argsort()[::-1]] 164 sorted_rt = time_subset[eic_subset.argsort()[::-1]] 165 166 # Calculate the dispersity index 167 cum_sum = np.cumsum(sorted_eic) / np.sum(sorted_eic) 168 rt_summ = sorted_rt[np.where(cum_sum < 0.5)] 169 if len(rt_summ) > 1: 170 d = np.std(rt_summ) 171 self._dispersity_index = d # minutes 172 self._normalized_dispersity_index = d / ( 173 time_subset[-1] - time_subset[0] 174 ) # unitless (fraction of total window used) 175 elif len(rt_summ) == 1: 176 self._dispersity_index = 0 177 self._normalized_dispersity_index = 0
Calculate the dispersity index of the mass feature.
This function calculates the dispersity index of the mass feature and
stores the result in the _dispersity_index attribute. The dispersity index is calculated as the standard
deviation of the retention times that account for 50% of the cummulative intensity, starting from the most
intense point, as described in [1]. Note that this calculation is done within the integration bounds with
a pad according to the window factor, where the window factor is parameterized and encapsulated in the
parent LCMS object (or, if not available, defaults to 2.0 minutes before and after the apex
Returns
- None, stores the result in the
_dispersity_indexattribute of the class and the_normalized_dispersity_indexattribute, - which is the dispersity index normalized to the total time window used for the calculation (unitless, fraction of total window).
Raises
- ValueError: If the EIC data are not available.
References
1) Boiteau, Rene M., et al. "Relating Molecular Properties to the Persistence of Marine Dissolved Organic Matter with Liquid Chromatography–Ultrahigh-Resolution Mass Spectrometry." Environmental Science & Technology 58.7 (2024): 3267-3277.
def
calc_fraction_height_width(self, fraction: float):
179 def calc_fraction_height_width(self, fraction: float): 180 """ 181 Calculate the height width of the mass feature at a specfic fraction of the maximum intensity. 182 183 This function returns a tuple with the minimum and maximum half-height width based on scan resolution. 184 185 Parameters 186 ---------- 187 fraction : float 188 The fraction of the maximum intensity to calculate the height width. 189 For example, 0.5 will calculate the half-height width. 190 191 Returns 192 ------- 193 Tuple[float, float, bool] 194 The minimum and maximum half-height width based on scan resolution (in minutes), and a boolean indicating if the width was estimated. 195 """ 196 197 # Pull out the EIC data 198 eic = self._eic_data.eic_smoothed 199 200 # Find the indices of the maximum intensity on either side 201 max_index = np.where(self._eic_data.scans == self.apex_scan)[0][0] 202 left_index = max_index 203 right_index = max_index 204 while eic[left_index] > eic[max_index] * fraction and left_index > 0: 205 left_index -= 1 206 while ( 207 eic[right_index] > eic[max_index] * fraction and right_index < len(eic) - 1 208 ): 209 right_index += 1 210 211 # Get the retention times of the indexes just below the half height 212 left_rt = self._eic_data.time[left_index] 213 right_rt = self._eic_data.time[right_index] 214 215 # If left_rt and right_rt are outside the bounds of the integration, set them to the bounds and set estimated to True 216 estimated = False 217 if left_rt < self.eic_rt_list[0]: 218 left_rt = self.eic_rt_list[0] 219 left_index = np.where(self._eic_data.scans == self._eic_data.apexes[0][0])[ 220 0 221 ][0] 222 estimated = True 223 if right_rt > self.eic_rt_list[-1]: 224 right_rt = self.eic_rt_list[-1] 225 right_index = np.where( 226 self._eic_data.scans == self._eic_data.apexes[0][-1] 227 )[0][0] 228 estimated = True 229 half_height_width_max = right_rt - left_rt 230 231 # Get the retention times of the indexes just above the half height 232 left_rt = self._eic_data.time[left_index + 1] 233 right_rt = self._eic_data.time[right_index - 1] 234 half_height_width_min = right_rt - left_rt 235 236 return half_height_width_min, half_height_width_max, estimated
Calculate the height width of the mass feature at a specfic fraction of the maximum intensity.
This function returns a tuple with the minimum and maximum half-height width based on scan resolution.
Parameters
- fraction (float): The fraction of the maximum intensity to calculate the height width. For example, 0.5 will calculate the half-height width.
Returns
- Tuple[float, float, bool]: The minimum and maximum half-height width based on scan resolution (in minutes), and a boolean indicating if the width was estimated.
def
calc_half_height_width(self, accept_estimated: bool = False):
238 def calc_half_height_width(self, accept_estimated: bool = False): 239 """ 240 Calculate the half-height width of the mass feature. 241 242 This function calculates the half-height width of the mass feature and 243 stores the result in the `_half_height_width` attribute 244 245 Returns 246 ------- 247 None, stores the result in the `_half_height_width` attribute of the class. 248 """ 249 min_, max_, estimated = self.calc_fraction_height_width(0.5) 250 if not estimated or accept_estimated: 251 self._half_height_width = np.array([min_, max_])
Calculate the half-height width of the mass feature.
This function calculates the half-height width of the mass feature and
stores the result in the _half_height_width attribute
Returns
- None, stores the result in the
_half_height_widthattribute of the class.
def
calc_tailing_factor(self, accept_estimated: bool = False):
253 def calc_tailing_factor(self, accept_estimated: bool = False): 254 """ 255 Calculate the peak asymmetry of the mass feature. 256 257 This function calculates the peak asymmetry of the mass feature and 258 stores the result in the `_tailing_factor` attribute. 259 Calculations completed at 5% of the peak height in accordance with the USP tailing factor calculation. 260 261 Returns 262 ------- 263 None, stores the result in the `_tailing_factor` attribute of the class. 264 265 References 266 ---------- 267 1) JIS K0124:2011 General rules for high performance liquid chromatography 268 2) JIS K0214:2013 Technical terms for analytical chemistry 269 """ 270 # First calculate the width of the peak at 5% of the peak height 271 width_min, width_max, estimated = self.calc_fraction_height_width(0.05) 272 273 if not estimated or accept_estimated: 274 # Next calculate the width of the peak at 95% of the peak height 275 eic = self._eic_data.eic_smoothed 276 max_index = np.where(self._eic_data.scans == self.apex_scan)[0][0] 277 left_index = max_index 278 while eic[left_index] > eic[max_index] * 0.05 and left_index > 0: 279 left_index -= 1 280 281 left_half_time_min = ( 282 self._eic_data.time[max_index] - self._eic_data.time[left_index] 283 ) 284 left_half_time_max = ( 285 self._eic_data.time[max_index] - self._eic_data.time[left_index + 1] 286 ) 287 288 tailing_factor = np.mean([width_min, width_max]) / ( 289 2 * np.mean([left_half_time_min, left_half_time_max]) 290 ) 291 292 self._tailing_factor = tailing_factor
Calculate the peak asymmetry of the mass feature.
This function calculates the peak asymmetry of the mass feature and
stores the result in the _tailing_factor attribute.
Calculations completed at 5% of the peak height in accordance with the USP tailing factor calculation.
Returns
- None, stores the result in the
_tailing_factorattribute of the class.
References
1) JIS K0124:2011 General rules for high performance liquid chromatography 2) JIS K0214:2013 Technical terms for analytical chemistry
def
calc_gaussian_similarity(self):
294 def calc_gaussian_similarity(self): 295 """ 296 Calculate the Gaussian similarity score of the mass feature. 297 298 This function fits a Gaussian curve to the EIC data and evaluates 299 the goodness of fit using R-squared. Note that this only uses data within 300 the set integration bounds of the mass feature. A score close to 1 indicates 301 the peak closely resembles an ideal Gaussian shape. 302 303 Returns 304 ------- 305 None, stores the result in the `_gaussian_similarity` attribute of the class. 306 307 Raises 308 ------ 309 ValueError 310 If the EIC data are not available. 311 """ 312 # Check if the EIC data is available 313 if self.eic_list is None: 314 raise ValueError( 315 "EIC data are not available. Please add the EIC data first." 316 ) 317 318 # Get EIC data within integration bounds 319 time_data = np.array(self.eic_rt_list) 320 intensity_data = np.array(self.eic_list) 321 322 if len(time_data) < 4: # Need minimum points for meaningful fit 323 self._gaussian_similarity = np.nan 324 return 325 326 # Check for valid intensity data 327 max_intensity = np.max(intensity_data) 328 if max_intensity == 0: 329 self._gaussian_similarity = np.nan 330 return 331 332 try: 333 # Define Gaussian function 334 def gaussian(x, amplitude, mean, stddev, baseline): 335 return ( 336 amplitude * np.exp(-((x - mean) ** 2) / (2 * stddev**2)) + baseline 337 ) 338 339 # Initial parameter estimates 340 amplitude_init = max_intensity 341 mean_init = time_data[np.argmax(intensity_data)] 342 stddev_init = (time_data[-1] - time_data[0]) / 6 # Rough estimate 343 baseline_init = np.min(intensity_data) 344 345 # Fit Gaussian curve 346 popt, _ = curve_fit( 347 gaussian, 348 time_data, 349 intensity_data, 350 p0=[amplitude_init, mean_init, stddev_init, baseline_init], 351 maxfev=1000, 352 bounds=( 353 [0, time_data[0], 0, 0], # Lower bounds 354 [np.inf, time_data[-1], np.inf, max_intensity], # Upper bounds 355 ), 356 ) 357 358 # Calculate fitted values 359 fitted_intensities = gaussian(time_data, *popt) 360 361 # Calculate R-squared (coefficient of determination) 362 ss_res = np.sum((intensity_data - fitted_intensities) ** 2) 363 ss_tot = np.sum((intensity_data - np.mean(intensity_data)) ** 2) 364 365 if ss_tot == 0: 366 self._gaussian_similarity = np.nan 367 else: 368 r_squared = 1 - (ss_res / ss_tot) 369 # R² should be between 0 and 1 for reasonable fits 370 # If negative, the model is worse than the mean - treat as non-computable 371 self._gaussian_similarity = r_squared if r_squared >= 0 else np.nan 372 373 except (RuntimeError, ValueError, TypeError): 374 # Fitting failed, assign NaN 375 self._gaussian_similarity = np.nan
Calculate the Gaussian similarity score of the mass feature.
This function fits a Gaussian curve to the EIC data and evaluates the goodness of fit using R-squared. Note that this only uses data within the set integration bounds of the mass feature. A score close to 1 indicates the peak closely resembles an ideal Gaussian shape.
Returns
- None, stores the result in the
_gaussian_similarityattribute of the class.
Raises
- ValueError: If the EIC data are not available.
def
calc_noise_score(self):
377 def calc_noise_score(self): 378 """ 379 Calculate the noise score of the mass feature separately for left and right sides. 380 381 This function estimates the signal-to-noise ratio by comparing the peak 382 intensity to the baseline noise level in surrounding regions. It calculates 383 separate scores for the left and right sides of the peak, which are stored as a tuple 384 in the `_noise_score` attribute. The noise estimation windows are encapsulated in the 385 parent LCMS object (or, if not available, defaults to twice the peak width on each side). 386 387 388 Returns 389 ------- 390 None, stores the result in the `_noise_score` attribute as a tuple (left_score, right_score). 391 392 Raises 393 ------ 394 ValueError 395 If the EIC data are not available. 396 """ 397 # Check if the EIC data is available 398 if self.eic_list is None: 399 raise ValueError( 400 "EIC data are not available. Please add the EIC data first." 401 ) 402 403 # Check if LCMSMassFeature has a parent LCMS object with a window factor 404 if hasattr(self, "mass_spectrum_obj"): 405 noise_window_factor = ( 406 self.mass_spectrum_obj.parameters.lc_ms.noise_window_factor 407 ) 408 else: 409 noise_window_factor = 2.0 # times the peak width 410 411 # Get full EIC data (not just integration bounds) 412 full_time = self._eic_data.time 413 full_eic = self._eic_data.eic 414 415 # Get peak information 416 apex_rt = self.retention_time 417 peak_intensity = np.max(self.eic_list) 418 419 # Retrieve width from integration bounds 420 peak_width = self.eic_rt_list[-1] - self.eic_rt_list[0] 421 422 # Define noise estimation windows 423 noise_window_size = peak_width * noise_window_factor # in minutes 424 left_noise_start = apex_rt - peak_width - noise_window_size 425 left_noise_end = apex_rt - peak_width 426 right_noise_start = apex_rt + peak_width 427 right_noise_end = apex_rt + peak_width + noise_window_size 428 429 # Extract noise regions 430 left_noise_mask = (full_time >= left_noise_start) & ( 431 full_time <= left_noise_end 432 ) 433 right_noise_mask = (full_time >= right_noise_start) & ( 434 full_time <= right_noise_end 435 ) 436 437 left_noise = full_eic[left_noise_mask] 438 right_noise = full_eic[right_noise_mask] 439 440 # Calculate left noise score 441 if len(left_noise) == 0: 442 left_score = np.nan 443 else: 444 left_baseline = np.median(left_noise) 445 left_noise_std = np.std(left_noise) 446 447 if left_noise_std == 0: 448 if peak_intensity > left_baseline: 449 left_score = 1.0 450 else: 451 left_score = np.nan 452 else: 453 left_signal = peak_intensity - left_baseline 454 if left_signal <= 0: 455 left_score = 0.0 456 else: 457 left_snr = left_signal / left_noise_std 458 left_score = min(1.0, left_snr / (left_snr + 10.0)) 459 460 # Calculate right noise score 461 if len(right_noise) == 0: 462 right_score = np.nan 463 else: 464 right_baseline = np.median(right_noise) 465 right_noise_std = np.std(right_noise) 466 467 if right_noise_std == 0: 468 if peak_intensity > right_baseline: 469 right_score = 1.0 470 else: 471 right_score = np.nan 472 else: 473 right_signal = peak_intensity - right_baseline 474 if right_signal <= 0: 475 right_score = 0.0 476 else: 477 right_snr = right_signal / right_noise_std 478 right_score = min(1.0, right_snr / (right_snr + 10.0)) 479 480 # Store as tuple 481 self._noise_score = (left_score, right_score)
Calculate the noise score of the mass feature separately for left and right sides.
This function estimates the signal-to-noise ratio by comparing the peak
intensity to the baseline noise level in surrounding regions. It calculates
separate scores for the left and right sides of the peak, which are stored as a tuple
in the _noise_score attribute. The noise estimation windows are encapsulated in the
parent LCMS object (or, if not available, defaults to twice the peak width on each side).
Returns
- None, stores the result in the
_noise_scoreattribute as a tuple (left_score, right_score).
Raises
- ValueError: If the EIC data are not available.