corems.mass_spectrum.calc.NoiseCalc
1import warnings 2from typing import Tuple 3 4from numpy import average, histogram, hstack, inf, isnan, log10, median, nan, std, where 5 6from corems import chunks 7 8# from matplotlib import pyplot 9__author__ = "Yuri E. Corilo" 10__date__ = "Jun 27, 2019" 11 12 13class NoiseThresholdCalc: 14 """Class for noise threshold calculation. 15 16 Parameters 17 ---------- 18 mass_spectrum : MassSpectrum 19 The mass spectrum object. 20 settings : MSParameters 21 The mass spectrum parameters object. 22 is_centroid : bool 23 Flag indicating whether the mass spectrum is centroid or profile. 24 baseline_noise : float 25 The baseline noise. 26 baseline_noise_std : float 27 The baseline noise standard deviation. 28 max_signal_to_noise : float 29 The maximum signal to noise. 30 max_abundance : float 31 The maximum abundance. 32 abundance : np.array 33 The abundance array. 34 abundance_profile : np.array 35 The abundance profile array. 36 mz_exp : np.array 37 The experimental m/z array. 38 mz_exp_profile : np.array 39 The experimental m/z profile array. 40 41 Attributes 42 ---------- 43 None 44 45 Methods 46 ------- 47 * get_noise_threshold(). Get the noise threshold. 48 * cut_mz_domain_noise(). Cut the m/z domain to the noise threshold regions. 49 * get_noise_average(ymincentroid). 50 Get the average noise and standard deviation. 51 * get_abundance_minima_centroid(abun_cut) 52 Get the abundance minima for centroid data. 53 * run_log_noise_threshold_calc(). 54 Run the log noise threshold calculation. 55 * run_noise_threshold_calc(). 56 Run the noise threshold calculation. 57 """ 58 59 def get_noise_threshold(self) -> Tuple[Tuple[float, float], Tuple[float, float]]: 60 """Get the noise threshold. 61 62 Returns 63 ------- 64 Tuple[Tuple[float, float], Tuple[float, float]] 65 A tuple containing the m/z and abundance noise thresholds. 66 (min_mz, max_mz), (noise_threshold, noise_threshold) 67 """ 68 69 self._resolve_noise_threshold_method() 70 71 if self.is_centroid: 72 x = min(self.mz_exp), max((self.mz_exp)) 73 74 if self.settings.noise_threshold_method == "minima": 75 abundance_threshold = self.baseline_noise + ( 76 self.settings.noise_threshold_min_std * self.baseline_noise_std 77 ) 78 y = (abundance_threshold, abundance_threshold) 79 80 elif self.settings.noise_threshold_method == "signal_noise": 81 normalized_threshold = ( 82 self.max_abundance * self.settings.noise_threshold_min_s2n 83 ) / self.max_signal_to_noise 84 y = (normalized_threshold, normalized_threshold) 85 86 elif self.settings.noise_threshold_method == "relative_abundance": 87 normalized_threshold = ( 88 max(self.abundance) / 100 89 ) * self.settings.noise_threshold_min_relative_abundance 90 y = (normalized_threshold, normalized_threshold) 91 92 elif self.settings.noise_threshold_method == "absolute_abundance": 93 normalized_threshold = ( 94 self.abundance * self.settings.noise_threshold_absolute_abundance 95 ) 96 y = (normalized_threshold, normalized_threshold) 97 # log noise method not tested for centroid data 98 else: 99 raise Exception( 100 "%s method was not implemented, please refer to corems.mass_spectrum.calc.NoiseCalc Class" 101 % self.settings.noise_threshold_method 102 ) 103 104 return x, y 105 106 else: 107 if self.baseline_noise and self.baseline_noise_std: 108 x = (self.mz_exp_profile.min(), self.mz_exp_profile.max()) 109 y = (self.baseline_noise_std, self.baseline_noise_std) 110 111 if self.settings.noise_threshold_method == "minima": 112 # print(self.settings.noise_threshold_min_std) 113 abundance_threshold = self.baseline_noise + ( 114 self.settings.noise_threshold_min_std * self.baseline_noise_std 115 ) 116 117 y = (abundance_threshold, abundance_threshold) 118 119 elif self.settings.noise_threshold_method == "signal_noise": 120 max_sn = self.abundance_profile.max() / self.baseline_noise_std 121 122 normalized_threshold = ( 123 self.abundance_profile.max() 124 * self.settings.noise_threshold_min_s2n 125 ) / max_sn 126 y = (normalized_threshold, normalized_threshold) 127 128 elif self.settings.noise_threshold_method == "relative_abundance": 129 normalized_threshold = ( 130 self.abundance_profile.max() / 100 131 ) * self.settings.noise_threshold_min_relative_abundance 132 y = (normalized_threshold, normalized_threshold) 133 134 elif self.settings.noise_threshold_method == "absolute_abundance": 135 normalized_threshold = ( 136 self.settings.noise_threshold_absolute_abundance 137 ) 138 y = (normalized_threshold, normalized_threshold) 139 140 elif self.settings.noise_threshold_method == "log": 141 normalized_threshold = ( 142 self.settings.noise_threshold_log_nsigma 143 * self.baseline_noise_std 144 ) 145 y = (normalized_threshold, normalized_threshold) 146 147 else: 148 raise Exception( 149 "%s method was not implemented, \ 150 please refer to corems.mass_spectrum.calc.NoiseCalc Class" 151 % self.settings.noise_threshold_method 152 ) 153 154 return x, y 155 156 else: 157 warnings.warn( 158 "Noise Baseline and Noise std not specified,\ 159 defaulting to 0,0 run process_mass_spec() ?" 160 ) 161 return (0, 0), (0, 0) 162 163 def cut_mz_domain_noise(self): 164 """Cut the m/z domain to the noise threshold regions. 165 166 Returns 167 ------- 168 Tuple[np.array, np.array] 169 A tuple containing the m/z and abundance arrays of the truncated spectrum region. 170 """ 171 min_mz_whole_ms = self.mz_exp_profile.min() 172 max_mz_whole_ms = self.mz_exp_profile.max() 173 174 if self.settings.noise_threshold_method == "minima": 175 # this calculation is taking too long (about 2 seconds) 176 number_average_molecular_weight = self.weight_average_molecular_weight( 177 profile=True 178 ) 179 180 # +-200 is a guess for testing only, it needs adjustment for each type of analysis 181 # need to check min mz here or it will break 182 min_mz_noise = number_average_molecular_weight - 100 183 # need to check max mz here or it will break 184 max_mz_noise = number_average_molecular_weight + 100 185 186 else: 187 min_mz_noise = self.settings.noise_min_mz 188 max_mz_noise = self.settings.noise_max_mz 189 190 if min_mz_noise < min_mz_whole_ms: 191 min_mz_noise = min_mz_whole_ms 192 193 if max_mz_noise > max_mz_whole_ms: 194 max_mz_noise = max_mz_whole_ms 195 196 if min_mz_noise > max_mz_noise: 197 warnings.warn( 198 "Empty noise ROI: min_mz_noise is greater than max_mz_noise; returning empty arrays." 199 ) 200 return self.mz_exp_profile[0:0], self.abundance_profile[0:0] 201 202 # Inclusive ROI boundaries to preserve current threshold behavior. 203 mz_mask = (self.mz_exp_profile >= min_mz_noise) & ( 204 self.mz_exp_profile <= max_mz_noise 205 ) 206 indices = where(mz_mask)[0] 207 208 if indices.size == 0: 209 warnings.warn( 210 "Empty noise ROI: no m/z points found within inclusive range; returning empty arrays." 211 ) 212 return self.mz_exp_profile[0:0], self.abundance_profile[0:0] 213 214 start_index = indices.min() 215 end_index = indices.max() + 1 216 217 return ( 218 self.mz_exp_profile[start_index:end_index], 219 self.abundance_profile[start_index:end_index], 220 ) 221 222 def get_noise_average(self, ymincentroid): 223 """Get the average noise and standard deviation. 224 225 Parameters 226 ---------- 227 ymincentroid : np.array 228 The ymincentroid array. 229 230 Returns 231 ------- 232 Tuple[float, float] 233 A tuple containing the average noise and standard deviation. 234 235 """ 236 # assumes noise to be gaussian and estimate noise level by 237 # calculating the valley. 238 239 auto = True if self.settings.noise_threshold_method == "minima" else False 240 241 average_noise = median((ymincentroid)) * 2 if auto else median(ymincentroid) 242 243 s_deviation = ymincentroid.std() * 3 if auto else ymincentroid.std() 244 245 return average_noise, s_deviation 246 247 def get_abundance_minima_centroid(self, abun_cut): 248 """Get the abundance minima for centroid data. 249 250 Parameters 251 ---------- 252 abun_cut : np.array 253 The abundance cut array. 254 255 Returns 256 ------- 257 np.array 258 The abundance minima array. 259 """ 260 maximum = self.abundance_profile.max() 261 threshold_min = maximum * 1.00 262 263 y = -abun_cut 264 265 dy = y[1:] - y[:-1] 266 """replaces NaN for Infinity""" 267 indices_nan = where(isnan(y))[0] 268 269 if indices_nan.size: 270 y[indices_nan] = inf 271 dy[where(isnan(dy))[0]] = inf 272 273 indices = where((hstack((dy, 0)) < 0) & (hstack((0, dy)) > 0))[0] 274 275 if indices.size and threshold_min is not None: 276 indices = indices[abun_cut[indices] <= threshold_min] 277 278 return abun_cut[indices] 279 280 def run_log_noise_threshold_calc(self): 281 """Run the log noise threshold calculation. 282 283 284 Returns 285 ------- 286 Tuple[float, float] 287 A tuple containing the average noise and standard deviation. 288 289 Notes 290 -------- 291 Method for estimating the noise based on decimal log of all the data point 292 293 Idea is that you calculate a histogram of of the log10(abundance) values. 294 The maximum of the histogram == the standard deviation of the noise. 295 296 297 For aFT data it is a gaussian distribution of noise - not implemented here! 298 For mFT data it is a Rayleigh distribution, and the value is actually 10^(abu_max)*0.463. 299 300 301 See the publication cited above for the derivation of this. 302 303 References 304 -------- 305 1. dx.doi.org/10.1021/ac403278t | Anal. Chem. 2014, 86, 3308−3316 306 307 """ 308 309 if self.is_centroid: 310 raise Exception("log noise Not tested for centroid data") 311 else: 312 # cut the spectrum to ROI 313 mz_cut, abundance_cut = self.cut_mz_domain_noise() 314 # If there are 0 values, the log will fail 315 # But we may have negative values for aFT data, so we check if 0 exists 316 # Need to make a copy of the abundance cut values so we dont overwrite it.... 317 tmp_abundance = abundance_cut.copy() 318 if 0 in tmp_abundance: 319 tmp_abundance[tmp_abundance == 0] = nan 320 tmp_abundance = tmp_abundance[~isnan(tmp_abundance)] 321 # It seems there are edge cases of sparse but high S/N data where the wrong values may be determined. 322 # Hard to generalise - needs more investigation. 323 324 # calculate a histogram of the log10 of the abundance data 325 hist_values = histogram( 326 log10(tmp_abundance), bins=self.settings.noise_threshold_log_nsigma_bins 327 ) 328 # find the apex of this histogram 329 maxvalidx = where(hist_values[0] == max(hist_values[0])) 330 # get the value of this apex (note - still in log10 units) 331 log_sigma = hist_values[1][maxvalidx] 332 # If the histogram had more than one maximum frequency bin, we need to reduce that to one entry 333 if len(log_sigma) > 1: 334 log_sigma = average(log_sigma) 335 else: 336 log_sigma = log_sigma[0] 337 ## To do : check if aFT or mFT and adjust method 338 noise_mid = 10**log_sigma 339 noise_1std = ( 340 noise_mid * self.settings.noise_threshold_log_nsigma_corr_factor 341 ) # for mFT 0.463 342 return float(noise_mid), float(noise_1std) 343 344 def run_noise_threshold_calc(self): 345 """Runs noise threshold calculation (not log based method) 346 347 Returns 348 ------- 349 Tuple[float, float] 350 A tuple containing the average noise and standard deviation. 351 352 """ 353 if self.is_centroid: 354 # calculates noise_baseline and noise_std 355 # needed to run auto noise threshold mode 356 # it is not used for signal to noise nor 357 # relative abudance methods 358 abundances_chunks = chunks(self.abundance, 50) 359 each_min_abund = [min(x) for x in abundances_chunks] 360 361 return average(each_min_abund), std(each_min_abund) 362 363 else: 364 mz_cut, abundance_cut = self.cut_mz_domain_noise() 365 366 if self.settings.noise_threshold_method == "minima": 367 yminima = self.get_abundance_minima_centroid(abundance_cut) 368 369 return self.get_noise_average(yminima) 370 371 else: 372 # pyplot.show() 373 return self.get_noise_average(abundance_cut)
class
NoiseThresholdCalc:
14class NoiseThresholdCalc: 15 """Class for noise threshold calculation. 16 17 Parameters 18 ---------- 19 mass_spectrum : MassSpectrum 20 The mass spectrum object. 21 settings : MSParameters 22 The mass spectrum parameters object. 23 is_centroid : bool 24 Flag indicating whether the mass spectrum is centroid or profile. 25 baseline_noise : float 26 The baseline noise. 27 baseline_noise_std : float 28 The baseline noise standard deviation. 29 max_signal_to_noise : float 30 The maximum signal to noise. 31 max_abundance : float 32 The maximum abundance. 33 abundance : np.array 34 The abundance array. 35 abundance_profile : np.array 36 The abundance profile array. 37 mz_exp : np.array 38 The experimental m/z array. 39 mz_exp_profile : np.array 40 The experimental m/z profile array. 41 42 Attributes 43 ---------- 44 None 45 46 Methods 47 ------- 48 * get_noise_threshold(). Get the noise threshold. 49 * cut_mz_domain_noise(). Cut the m/z domain to the noise threshold regions. 50 * get_noise_average(ymincentroid). 51 Get the average noise and standard deviation. 52 * get_abundance_minima_centroid(abun_cut) 53 Get the abundance minima for centroid data. 54 * run_log_noise_threshold_calc(). 55 Run the log noise threshold calculation. 56 * run_noise_threshold_calc(). 57 Run the noise threshold calculation. 58 """ 59 60 def get_noise_threshold(self) -> Tuple[Tuple[float, float], Tuple[float, float]]: 61 """Get the noise threshold. 62 63 Returns 64 ------- 65 Tuple[Tuple[float, float], Tuple[float, float]] 66 A tuple containing the m/z and abundance noise thresholds. 67 (min_mz, max_mz), (noise_threshold, noise_threshold) 68 """ 69 70 self._resolve_noise_threshold_method() 71 72 if self.is_centroid: 73 x = min(self.mz_exp), max((self.mz_exp)) 74 75 if self.settings.noise_threshold_method == "minima": 76 abundance_threshold = self.baseline_noise + ( 77 self.settings.noise_threshold_min_std * self.baseline_noise_std 78 ) 79 y = (abundance_threshold, abundance_threshold) 80 81 elif self.settings.noise_threshold_method == "signal_noise": 82 normalized_threshold = ( 83 self.max_abundance * self.settings.noise_threshold_min_s2n 84 ) / self.max_signal_to_noise 85 y = (normalized_threshold, normalized_threshold) 86 87 elif self.settings.noise_threshold_method == "relative_abundance": 88 normalized_threshold = ( 89 max(self.abundance) / 100 90 ) * self.settings.noise_threshold_min_relative_abundance 91 y = (normalized_threshold, normalized_threshold) 92 93 elif self.settings.noise_threshold_method == "absolute_abundance": 94 normalized_threshold = ( 95 self.abundance * self.settings.noise_threshold_absolute_abundance 96 ) 97 y = (normalized_threshold, normalized_threshold) 98 # log noise method not tested for centroid data 99 else: 100 raise Exception( 101 "%s method was not implemented, please refer to corems.mass_spectrum.calc.NoiseCalc Class" 102 % self.settings.noise_threshold_method 103 ) 104 105 return x, y 106 107 else: 108 if self.baseline_noise and self.baseline_noise_std: 109 x = (self.mz_exp_profile.min(), self.mz_exp_profile.max()) 110 y = (self.baseline_noise_std, self.baseline_noise_std) 111 112 if self.settings.noise_threshold_method == "minima": 113 # print(self.settings.noise_threshold_min_std) 114 abundance_threshold = self.baseline_noise + ( 115 self.settings.noise_threshold_min_std * self.baseline_noise_std 116 ) 117 118 y = (abundance_threshold, abundance_threshold) 119 120 elif self.settings.noise_threshold_method == "signal_noise": 121 max_sn = self.abundance_profile.max() / self.baseline_noise_std 122 123 normalized_threshold = ( 124 self.abundance_profile.max() 125 * self.settings.noise_threshold_min_s2n 126 ) / max_sn 127 y = (normalized_threshold, normalized_threshold) 128 129 elif self.settings.noise_threshold_method == "relative_abundance": 130 normalized_threshold = ( 131 self.abundance_profile.max() / 100 132 ) * self.settings.noise_threshold_min_relative_abundance 133 y = (normalized_threshold, normalized_threshold) 134 135 elif self.settings.noise_threshold_method == "absolute_abundance": 136 normalized_threshold = ( 137 self.settings.noise_threshold_absolute_abundance 138 ) 139 y = (normalized_threshold, normalized_threshold) 140 141 elif self.settings.noise_threshold_method == "log": 142 normalized_threshold = ( 143 self.settings.noise_threshold_log_nsigma 144 * self.baseline_noise_std 145 ) 146 y = (normalized_threshold, normalized_threshold) 147 148 else: 149 raise Exception( 150 "%s method was not implemented, \ 151 please refer to corems.mass_spectrum.calc.NoiseCalc Class" 152 % self.settings.noise_threshold_method 153 ) 154 155 return x, y 156 157 else: 158 warnings.warn( 159 "Noise Baseline and Noise std not specified,\ 160 defaulting to 0,0 run process_mass_spec() ?" 161 ) 162 return (0, 0), (0, 0) 163 164 def cut_mz_domain_noise(self): 165 """Cut the m/z domain to the noise threshold regions. 166 167 Returns 168 ------- 169 Tuple[np.array, np.array] 170 A tuple containing the m/z and abundance arrays of the truncated spectrum region. 171 """ 172 min_mz_whole_ms = self.mz_exp_profile.min() 173 max_mz_whole_ms = self.mz_exp_profile.max() 174 175 if self.settings.noise_threshold_method == "minima": 176 # this calculation is taking too long (about 2 seconds) 177 number_average_molecular_weight = self.weight_average_molecular_weight( 178 profile=True 179 ) 180 181 # +-200 is a guess for testing only, it needs adjustment for each type of analysis 182 # need to check min mz here or it will break 183 min_mz_noise = number_average_molecular_weight - 100 184 # need to check max mz here or it will break 185 max_mz_noise = number_average_molecular_weight + 100 186 187 else: 188 min_mz_noise = self.settings.noise_min_mz 189 max_mz_noise = self.settings.noise_max_mz 190 191 if min_mz_noise < min_mz_whole_ms: 192 min_mz_noise = min_mz_whole_ms 193 194 if max_mz_noise > max_mz_whole_ms: 195 max_mz_noise = max_mz_whole_ms 196 197 if min_mz_noise > max_mz_noise: 198 warnings.warn( 199 "Empty noise ROI: min_mz_noise is greater than max_mz_noise; returning empty arrays." 200 ) 201 return self.mz_exp_profile[0:0], self.abundance_profile[0:0] 202 203 # Inclusive ROI boundaries to preserve current threshold behavior. 204 mz_mask = (self.mz_exp_profile >= min_mz_noise) & ( 205 self.mz_exp_profile <= max_mz_noise 206 ) 207 indices = where(mz_mask)[0] 208 209 if indices.size == 0: 210 warnings.warn( 211 "Empty noise ROI: no m/z points found within inclusive range; returning empty arrays." 212 ) 213 return self.mz_exp_profile[0:0], self.abundance_profile[0:0] 214 215 start_index = indices.min() 216 end_index = indices.max() + 1 217 218 return ( 219 self.mz_exp_profile[start_index:end_index], 220 self.abundance_profile[start_index:end_index], 221 ) 222 223 def get_noise_average(self, ymincentroid): 224 """Get the average noise and standard deviation. 225 226 Parameters 227 ---------- 228 ymincentroid : np.array 229 The ymincentroid array. 230 231 Returns 232 ------- 233 Tuple[float, float] 234 A tuple containing the average noise and standard deviation. 235 236 """ 237 # assumes noise to be gaussian and estimate noise level by 238 # calculating the valley. 239 240 auto = True if self.settings.noise_threshold_method == "minima" else False 241 242 average_noise = median((ymincentroid)) * 2 if auto else median(ymincentroid) 243 244 s_deviation = ymincentroid.std() * 3 if auto else ymincentroid.std() 245 246 return average_noise, s_deviation 247 248 def get_abundance_minima_centroid(self, abun_cut): 249 """Get the abundance minima for centroid data. 250 251 Parameters 252 ---------- 253 abun_cut : np.array 254 The abundance cut array. 255 256 Returns 257 ------- 258 np.array 259 The abundance minima array. 260 """ 261 maximum = self.abundance_profile.max() 262 threshold_min = maximum * 1.00 263 264 y = -abun_cut 265 266 dy = y[1:] - y[:-1] 267 """replaces NaN for Infinity""" 268 indices_nan = where(isnan(y))[0] 269 270 if indices_nan.size: 271 y[indices_nan] = inf 272 dy[where(isnan(dy))[0]] = inf 273 274 indices = where((hstack((dy, 0)) < 0) & (hstack((0, dy)) > 0))[0] 275 276 if indices.size and threshold_min is not None: 277 indices = indices[abun_cut[indices] <= threshold_min] 278 279 return abun_cut[indices] 280 281 def run_log_noise_threshold_calc(self): 282 """Run the log noise threshold calculation. 283 284 285 Returns 286 ------- 287 Tuple[float, float] 288 A tuple containing the average noise and standard deviation. 289 290 Notes 291 -------- 292 Method for estimating the noise based on decimal log of all the data point 293 294 Idea is that you calculate a histogram of of the log10(abundance) values. 295 The maximum of the histogram == the standard deviation of the noise. 296 297 298 For aFT data it is a gaussian distribution of noise - not implemented here! 299 For mFT data it is a Rayleigh distribution, and the value is actually 10^(abu_max)*0.463. 300 301 302 See the publication cited above for the derivation of this. 303 304 References 305 -------- 306 1. dx.doi.org/10.1021/ac403278t | Anal. Chem. 2014, 86, 3308−3316 307 308 """ 309 310 if self.is_centroid: 311 raise Exception("log noise Not tested for centroid data") 312 else: 313 # cut the spectrum to ROI 314 mz_cut, abundance_cut = self.cut_mz_domain_noise() 315 # If there are 0 values, the log will fail 316 # But we may have negative values for aFT data, so we check if 0 exists 317 # Need to make a copy of the abundance cut values so we dont overwrite it.... 318 tmp_abundance = abundance_cut.copy() 319 if 0 in tmp_abundance: 320 tmp_abundance[tmp_abundance == 0] = nan 321 tmp_abundance = tmp_abundance[~isnan(tmp_abundance)] 322 # It seems there are edge cases of sparse but high S/N data where the wrong values may be determined. 323 # Hard to generalise - needs more investigation. 324 325 # calculate a histogram of the log10 of the abundance data 326 hist_values = histogram( 327 log10(tmp_abundance), bins=self.settings.noise_threshold_log_nsigma_bins 328 ) 329 # find the apex of this histogram 330 maxvalidx = where(hist_values[0] == max(hist_values[0])) 331 # get the value of this apex (note - still in log10 units) 332 log_sigma = hist_values[1][maxvalidx] 333 # If the histogram had more than one maximum frequency bin, we need to reduce that to one entry 334 if len(log_sigma) > 1: 335 log_sigma = average(log_sigma) 336 else: 337 log_sigma = log_sigma[0] 338 ## To do : check if aFT or mFT and adjust method 339 noise_mid = 10**log_sigma 340 noise_1std = ( 341 noise_mid * self.settings.noise_threshold_log_nsigma_corr_factor 342 ) # for mFT 0.463 343 return float(noise_mid), float(noise_1std) 344 345 def run_noise_threshold_calc(self): 346 """Runs noise threshold calculation (not log based method) 347 348 Returns 349 ------- 350 Tuple[float, float] 351 A tuple containing the average noise and standard deviation. 352 353 """ 354 if self.is_centroid: 355 # calculates noise_baseline and noise_std 356 # needed to run auto noise threshold mode 357 # it is not used for signal to noise nor 358 # relative abudance methods 359 abundances_chunks = chunks(self.abundance, 50) 360 each_min_abund = [min(x) for x in abundances_chunks] 361 362 return average(each_min_abund), std(each_min_abund) 363 364 else: 365 mz_cut, abundance_cut = self.cut_mz_domain_noise() 366 367 if self.settings.noise_threshold_method == "minima": 368 yminima = self.get_abundance_minima_centroid(abundance_cut) 369 370 return self.get_noise_average(yminima) 371 372 else: 373 # pyplot.show() 374 return self.get_noise_average(abundance_cut)
Class for noise threshold calculation.
Parameters
- mass_spectrum (MassSpectrum): The mass spectrum object.
- settings (MSParameters): The mass spectrum parameters object.
- is_centroid (bool): Flag indicating whether the mass spectrum is centroid or profile.
- baseline_noise (float): The baseline noise.
- baseline_noise_std (float): The baseline noise standard deviation.
- max_signal_to_noise (float): The maximum signal to noise.
- max_abundance (float): The maximum abundance.
- abundance (np.array): The abundance array.
- abundance_profile (np.array): The abundance profile array.
- mz_exp (np.array): The experimental m/z array.
- mz_exp_profile (np.array): The experimental m/z profile array.
Attributes
- None
Methods
- get_noise_threshold(). Get the noise threshold.
- cut_mz_domain_noise(). Cut the m/z domain to the noise threshold regions.
- get_noise_average(ymincentroid). Get the average noise and standard deviation.
- get_abundance_minima_centroid(abun_cut) Get the abundance minima for centroid data.
- run_log_noise_threshold_calc(). Run the log noise threshold calculation.
- run_noise_threshold_calc(). Run the noise threshold calculation.
def
get_noise_threshold(self) -> Tuple[Tuple[float, float], Tuple[float, float]]:
60 def get_noise_threshold(self) -> Tuple[Tuple[float, float], Tuple[float, float]]: 61 """Get the noise threshold. 62 63 Returns 64 ------- 65 Tuple[Tuple[float, float], Tuple[float, float]] 66 A tuple containing the m/z and abundance noise thresholds. 67 (min_mz, max_mz), (noise_threshold, noise_threshold) 68 """ 69 70 self._resolve_noise_threshold_method() 71 72 if self.is_centroid: 73 x = min(self.mz_exp), max((self.mz_exp)) 74 75 if self.settings.noise_threshold_method == "minima": 76 abundance_threshold = self.baseline_noise + ( 77 self.settings.noise_threshold_min_std * self.baseline_noise_std 78 ) 79 y = (abundance_threshold, abundance_threshold) 80 81 elif self.settings.noise_threshold_method == "signal_noise": 82 normalized_threshold = ( 83 self.max_abundance * self.settings.noise_threshold_min_s2n 84 ) / self.max_signal_to_noise 85 y = (normalized_threshold, normalized_threshold) 86 87 elif self.settings.noise_threshold_method == "relative_abundance": 88 normalized_threshold = ( 89 max(self.abundance) / 100 90 ) * self.settings.noise_threshold_min_relative_abundance 91 y = (normalized_threshold, normalized_threshold) 92 93 elif self.settings.noise_threshold_method == "absolute_abundance": 94 normalized_threshold = ( 95 self.abundance * self.settings.noise_threshold_absolute_abundance 96 ) 97 y = (normalized_threshold, normalized_threshold) 98 # log noise method not tested for centroid data 99 else: 100 raise Exception( 101 "%s method was not implemented, please refer to corems.mass_spectrum.calc.NoiseCalc Class" 102 % self.settings.noise_threshold_method 103 ) 104 105 return x, y 106 107 else: 108 if self.baseline_noise and self.baseline_noise_std: 109 x = (self.mz_exp_profile.min(), self.mz_exp_profile.max()) 110 y = (self.baseline_noise_std, self.baseline_noise_std) 111 112 if self.settings.noise_threshold_method == "minima": 113 # print(self.settings.noise_threshold_min_std) 114 abundance_threshold = self.baseline_noise + ( 115 self.settings.noise_threshold_min_std * self.baseline_noise_std 116 ) 117 118 y = (abundance_threshold, abundance_threshold) 119 120 elif self.settings.noise_threshold_method == "signal_noise": 121 max_sn = self.abundance_profile.max() / self.baseline_noise_std 122 123 normalized_threshold = ( 124 self.abundance_profile.max() 125 * self.settings.noise_threshold_min_s2n 126 ) / max_sn 127 y = (normalized_threshold, normalized_threshold) 128 129 elif self.settings.noise_threshold_method == "relative_abundance": 130 normalized_threshold = ( 131 self.abundance_profile.max() / 100 132 ) * self.settings.noise_threshold_min_relative_abundance 133 y = (normalized_threshold, normalized_threshold) 134 135 elif self.settings.noise_threshold_method == "absolute_abundance": 136 normalized_threshold = ( 137 self.settings.noise_threshold_absolute_abundance 138 ) 139 y = (normalized_threshold, normalized_threshold) 140 141 elif self.settings.noise_threshold_method == "log": 142 normalized_threshold = ( 143 self.settings.noise_threshold_log_nsigma 144 * self.baseline_noise_std 145 ) 146 y = (normalized_threshold, normalized_threshold) 147 148 else: 149 raise Exception( 150 "%s method was not implemented, \ 151 please refer to corems.mass_spectrum.calc.NoiseCalc Class" 152 % self.settings.noise_threshold_method 153 ) 154 155 return x, y 156 157 else: 158 warnings.warn( 159 "Noise Baseline and Noise std not specified,\ 160 defaulting to 0,0 run process_mass_spec() ?" 161 ) 162 return (0, 0), (0, 0)
Get the noise threshold.
Returns
- Tuple[Tuple[float, float], Tuple[float, float]]: A tuple containing the m/z and abundance noise thresholds. (min_mz, max_mz), (noise_threshold, noise_threshold)
def
cut_mz_domain_noise(self):
164 def cut_mz_domain_noise(self): 165 """Cut the m/z domain to the noise threshold regions. 166 167 Returns 168 ------- 169 Tuple[np.array, np.array] 170 A tuple containing the m/z and abundance arrays of the truncated spectrum region. 171 """ 172 min_mz_whole_ms = self.mz_exp_profile.min() 173 max_mz_whole_ms = self.mz_exp_profile.max() 174 175 if self.settings.noise_threshold_method == "minima": 176 # this calculation is taking too long (about 2 seconds) 177 number_average_molecular_weight = self.weight_average_molecular_weight( 178 profile=True 179 ) 180 181 # +-200 is a guess for testing only, it needs adjustment for each type of analysis 182 # need to check min mz here or it will break 183 min_mz_noise = number_average_molecular_weight - 100 184 # need to check max mz here or it will break 185 max_mz_noise = number_average_molecular_weight + 100 186 187 else: 188 min_mz_noise = self.settings.noise_min_mz 189 max_mz_noise = self.settings.noise_max_mz 190 191 if min_mz_noise < min_mz_whole_ms: 192 min_mz_noise = min_mz_whole_ms 193 194 if max_mz_noise > max_mz_whole_ms: 195 max_mz_noise = max_mz_whole_ms 196 197 if min_mz_noise > max_mz_noise: 198 warnings.warn( 199 "Empty noise ROI: min_mz_noise is greater than max_mz_noise; returning empty arrays." 200 ) 201 return self.mz_exp_profile[0:0], self.abundance_profile[0:0] 202 203 # Inclusive ROI boundaries to preserve current threshold behavior. 204 mz_mask = (self.mz_exp_profile >= min_mz_noise) & ( 205 self.mz_exp_profile <= max_mz_noise 206 ) 207 indices = where(mz_mask)[0] 208 209 if indices.size == 0: 210 warnings.warn( 211 "Empty noise ROI: no m/z points found within inclusive range; returning empty arrays." 212 ) 213 return self.mz_exp_profile[0:0], self.abundance_profile[0:0] 214 215 start_index = indices.min() 216 end_index = indices.max() + 1 217 218 return ( 219 self.mz_exp_profile[start_index:end_index], 220 self.abundance_profile[start_index:end_index], 221 )
Cut the m/z domain to the noise threshold regions.
Returns
- Tuple[np.array, np.array]: A tuple containing the m/z and abundance arrays of the truncated spectrum region.
def
get_noise_average(self, ymincentroid):
223 def get_noise_average(self, ymincentroid): 224 """Get the average noise and standard deviation. 225 226 Parameters 227 ---------- 228 ymincentroid : np.array 229 The ymincentroid array. 230 231 Returns 232 ------- 233 Tuple[float, float] 234 A tuple containing the average noise and standard deviation. 235 236 """ 237 # assumes noise to be gaussian and estimate noise level by 238 # calculating the valley. 239 240 auto = True if self.settings.noise_threshold_method == "minima" else False 241 242 average_noise = median((ymincentroid)) * 2 if auto else median(ymincentroid) 243 244 s_deviation = ymincentroid.std() * 3 if auto else ymincentroid.std() 245 246 return average_noise, s_deviation
Get the average noise and standard deviation.
Parameters
- ymincentroid (np.array): The ymincentroid array.
Returns
- Tuple[float, float]: A tuple containing the average noise and standard deviation.
def
get_abundance_minima_centroid(self, abun_cut):
248 def get_abundance_minima_centroid(self, abun_cut): 249 """Get the abundance minima for centroid data. 250 251 Parameters 252 ---------- 253 abun_cut : np.array 254 The abundance cut array. 255 256 Returns 257 ------- 258 np.array 259 The abundance minima array. 260 """ 261 maximum = self.abundance_profile.max() 262 threshold_min = maximum * 1.00 263 264 y = -abun_cut 265 266 dy = y[1:] - y[:-1] 267 """replaces NaN for Infinity""" 268 indices_nan = where(isnan(y))[0] 269 270 if indices_nan.size: 271 y[indices_nan] = inf 272 dy[where(isnan(dy))[0]] = inf 273 274 indices = where((hstack((dy, 0)) < 0) & (hstack((0, dy)) > 0))[0] 275 276 if indices.size and threshold_min is not None: 277 indices = indices[abun_cut[indices] <= threshold_min] 278 279 return abun_cut[indices]
Get the abundance minima for centroid data.
Parameters
- abun_cut (np.array): The abundance cut array.
Returns
- np.array: The abundance minima array.
def
run_log_noise_threshold_calc(self):
281 def run_log_noise_threshold_calc(self): 282 """Run the log noise threshold calculation. 283 284 285 Returns 286 ------- 287 Tuple[float, float] 288 A tuple containing the average noise and standard deviation. 289 290 Notes 291 -------- 292 Method for estimating the noise based on decimal log of all the data point 293 294 Idea is that you calculate a histogram of of the log10(abundance) values. 295 The maximum of the histogram == the standard deviation of the noise. 296 297 298 For aFT data it is a gaussian distribution of noise - not implemented here! 299 For mFT data it is a Rayleigh distribution, and the value is actually 10^(abu_max)*0.463. 300 301 302 See the publication cited above for the derivation of this. 303 304 References 305 -------- 306 1. dx.doi.org/10.1021/ac403278t | Anal. Chem. 2014, 86, 3308−3316 307 308 """ 309 310 if self.is_centroid: 311 raise Exception("log noise Not tested for centroid data") 312 else: 313 # cut the spectrum to ROI 314 mz_cut, abundance_cut = self.cut_mz_domain_noise() 315 # If there are 0 values, the log will fail 316 # But we may have negative values for aFT data, so we check if 0 exists 317 # Need to make a copy of the abundance cut values so we dont overwrite it.... 318 tmp_abundance = abundance_cut.copy() 319 if 0 in tmp_abundance: 320 tmp_abundance[tmp_abundance == 0] = nan 321 tmp_abundance = tmp_abundance[~isnan(tmp_abundance)] 322 # It seems there are edge cases of sparse but high S/N data where the wrong values may be determined. 323 # Hard to generalise - needs more investigation. 324 325 # calculate a histogram of the log10 of the abundance data 326 hist_values = histogram( 327 log10(tmp_abundance), bins=self.settings.noise_threshold_log_nsigma_bins 328 ) 329 # find the apex of this histogram 330 maxvalidx = where(hist_values[0] == max(hist_values[0])) 331 # get the value of this apex (note - still in log10 units) 332 log_sigma = hist_values[1][maxvalidx] 333 # If the histogram had more than one maximum frequency bin, we need to reduce that to one entry 334 if len(log_sigma) > 1: 335 log_sigma = average(log_sigma) 336 else: 337 log_sigma = log_sigma[0] 338 ## To do : check if aFT or mFT and adjust method 339 noise_mid = 10**log_sigma 340 noise_1std = ( 341 noise_mid * self.settings.noise_threshold_log_nsigma_corr_factor 342 ) # for mFT 0.463 343 return float(noise_mid), float(noise_1std)
Run the log noise threshold calculation.
Returns
- Tuple[float, float]: A tuple containing the average noise and standard deviation.
Notes
Method for estimating the noise based on decimal log of all the data point
Idea is that you calculate a histogram of of the log10(abundance) values. The maximum of the histogram == the standard deviation of the noise.
For aFT data it is a gaussian distribution of noise - not implemented here! For mFT data it is a Rayleigh distribution, and the value is actually 10^(abu_max)*0.463.
See the publication cited above for the derivation of this.
References
- dx.doi.org/10.1021/ac403278t | Anal. Chem. 2014, 86, 3308−3316
def
run_noise_threshold_calc(self):
345 def run_noise_threshold_calc(self): 346 """Runs noise threshold calculation (not log based method) 347 348 Returns 349 ------- 350 Tuple[float, float] 351 A tuple containing the average noise and standard deviation. 352 353 """ 354 if self.is_centroid: 355 # calculates noise_baseline and noise_std 356 # needed to run auto noise threshold mode 357 # it is not used for signal to noise nor 358 # relative abudance methods 359 abundances_chunks = chunks(self.abundance, 50) 360 each_min_abund = [min(x) for x in abundances_chunks] 361 362 return average(each_min_abund), std(each_min_abund) 363 364 else: 365 mz_cut, abundance_cut = self.cut_mz_domain_noise() 366 367 if self.settings.noise_threshold_method == "minima": 368 yminima = self.get_abundance_minima_centroid(abundance_cut) 369 370 return self.get_noise_average(yminima) 371 372 else: 373 # pyplot.show() 374 return self.get_noise_average(abundance_cut)
Runs noise threshold calculation (not log based method)
Returns
- Tuple[float, float]: A tuple containing the average noise and standard deviation.