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