corems.mass_spectrum.calc.MassErrorPrediction
1__author__ = "Yuri E. Corilo" 2__date__ = "03/31/2020" 3 4from threading import Thread 5from pandas import DataFrame 6from numpy import hstack, inf, isnan, where, array 7from tqdm import tqdm 8 9 10class MassErrorPrediction(Thread): 11 """Class for mass error prediction. 12 13 Parameters 14 ---------- 15 mass_spectrum : list 16 List of mass spectrum objects. 17 mz_overlay : int, optional 18 The mz overlay value for peak simulation. Default is 10. 19 rp_increments : int, optional 20 The resolving power increments for peak simulation. Default is 10000. 21 base_line_target : float, optional 22 The target value for the baseline resolution. Default is 0.01. 23 max_interation : int, optional 24 The maximum number of iterations for peak simulation. Default is 1000. 25 interpolation : str, optional 26 The interpolation method for missing data. Default is 'linear'. 27 28 Attributes 29 ---------- 30 mass_spectrum_obj : list 31 List of mass spectrum objects. 32 mz_overlay : int 33 The mz overlay value for peak simulation. 34 rp_increments : int 35 The resolving power increments for peak simulation. 36 base_line_target : float 37 The target value for the baseline resolution. 38 max_interation : int 39 The maximum number of iterations for peak simulation. 40 df : DataFrame or None 41 The calculated error distribution dataframe. 42 interpolation : str 43 The interpolation method for missing data. 44 45 Methods 46 ------- 47 * run(). 48 Runs the mass error prediction calculation. 49 * get_results(). 50 Returns the calculated error distribution dataframe. 51 52 """ 53 54 def __init__( 55 self, 56 mass_spectrum, 57 mz_overlay=10, 58 rp_increments=10000, 59 base_line_target: float = 0.01, 60 max_interation=1000, 61 interpolation="linear", 62 ): 63 Thread.__init__(self) 64 65 self.mass_spectrum_obj = mass_spectrum 66 67 self.mz_overlay = mz_overlay 68 69 self.rp_increments = rp_increments 70 71 self.base_line_target = base_line_target 72 73 self.max_interation = max_interation 74 75 self.df = None 76 77 self.interpolation = interpolation 78 79 def run(self): 80 """Runs the mass error prediction calculation.""" 81 self.df = self.calc_error_dist() 82 83 def get_results(self): 84 """Returns the calculated error distribution dataframe.""" 85 86 if not self.df: 87 self.run() 88 89 return self.df 90 91 def calc_error_dist(self): 92 """Calculate the error distribution.""" 93 verbose = self.mass_spectrum_obj.parameters.mass_spectrum.verbose_processing 94 results_list = [] 95 96 indexes_without_results = list(range(len(self.mass_spectrum_obj))) 97 # loop trough mass spectrum 98 for peak_obj_idx, peak_obj in enumerate( 99 tqdm(self.mass_spectrum_obj), disable=not verbose 100 ): 101 # access ms peaks triplets ( peak_obj_idx -1, peak_obj_idx, and peak_obj_idx + 1) 102 # check lower and upper boundaries to not excesses mass spectrum range 103 104 if peak_obj_idx != 0 and peak_obj_idx != len(self.mass_spectrum_obj) - 1: 105 # current peak_obj initialted in the loop expression 106 # geting the peak on the left (previous_peak_obj) and the one in the right position (next_peak_obj) 107 next_peak_obj = self.mass_spectrum_obj[peak_obj_idx + 1] 108 previous_peak_obj = self.mass_spectrum_obj[peak_obj_idx - 1] 109 110 # check mz range defined in max_mz variable and check if peaks have same nominal mz 111 # keeping same mz for better plotting representation only, remove it for production 112 if ( 113 peak_obj.nominal_mz_exp == next_peak_obj.nominal_mz_exp 114 and peak_obj.nominal_mz_exp == previous_peak_obj.nominal_mz_exp 115 ): 116 # simulate peak shape 117 sim_mz, sim_abun = peak_obj.gaussian(mz_overlay=self.mz_overlay) 118 # update_plot(sim_mz,sim_abun, 0.5) 119 120 # simulate peak shape 121 next_sim_mz, next_sim_abun = next_peak_obj.gaussian( 122 mz_overlay=self.mz_overlay 123 ) 124 # update_plot(next_sim_mz, next_sim_abun, 0.5) 125 126 # simulate peak shape 127 previous_sim_mz, previous_sim_abun = previous_peak_obj.gaussian( 128 mz_overlay=self.mz_overlay 129 ) 130 # update_plot(previous_sim_mz, previous_sim_abun, 0.5) 131 132 sim_mz_domain, summed_peaks_abun = self.sum_data( 133 ( 134 (previous_sim_mz, previous_sim_abun), 135 (sim_mz, sim_abun), 136 (next_sim_mz, next_sim_abun), 137 ) 138 ) 139 # update_plot(sim_mz_domain,summed_peaks_abun, 0.5) 140 141 # sum simulated abundances 142 # summed_peaks_abun = (sim_abun + next_sim_abun + previous_sim_abun) 143 144 # normalize abundances to 0-1 145 # summed_peaks_abun = summed_peaks_abun/(max(summed_peaks_abun)) 146 147 # find appexes location (mz) and magnitude 148 mz_centroid, abund_centroid = self.find_peak_apex( 149 sim_mz_domain, summed_peaks_abun 150 ) 151 152 # find valley location (mz_min_valley) and magnitude (abund_min_valley) 153 mz_min_valley, abund_min_valley = self.find_peak_valley( 154 sim_mz_domain, summed_peaks_abun 155 ) 156 157 # clear delta_rp (global implementation) and store choose resolving power increments 158 delta_rp = self.rp_increments 159 160 # used to limited number of iterations 161 i = 0 162 j = 0 163 164 # TODO: fit peak shape and decide best fit #gaussian, lorentz and voigt 165 # plot_triplets(mz_centroid,abund_centroid, mz_min_valley, abund_min_valley, sim_mz_domain, summed_peaks_abun ) 166 if len(mz_centroid) == 2: 167 while len(mz_centroid) < 3 and i <= self.max_interation: 168 previous_sim_mz, previous_sim_abun = ( 169 previous_peak_obj.gaussian( 170 delta_rp=delta_rp, mz_overlay=self.mz_overlay 171 ) 172 ) 173 174 sim_mz, sim_abun = peak_obj.gaussian( 175 delta_rp=delta_rp, mz_overlay=self.mz_overlay 176 ) 177 178 next_sim_mz, next_sim_abun = next_peak_obj.gaussian( 179 delta_rp=delta_rp, mz_overlay=self.mz_overlay 180 ) 181 182 sim_mz_domain, summed_peaks_abun = self.sum_data( 183 ( 184 (previous_sim_mz, previous_sim_abun), 185 (sim_mz, sim_abun), 186 (next_sim_mz, next_sim_abun), 187 ) 188 ) 189 190 # update_plot(sim_mz_domain, summed_peaks_abun, 0.01) 191 192 mz_centroid, abund_centroid = self.find_peak_apex( 193 sim_mz_domain, summed_peaks_abun 194 ) 195 196 delta_rp += self.rp_increments 197 198 i += 1 199 200 mz_min_valley, abund_min_valley = self.find_peak_valley( 201 sim_mz_domain, summed_peaks_abun 202 ) 203 204 if len(mz_centroid) == 3 and len(abund_min_valley) == 2: 205 # increase all three peak resolving power until both valley magnitude is bellow the defined target 206 # calculate peak shapes with the needed resolving power to have a baseline resolution for all peaks 207 # calculate mass difference (ppm) between original centroid and the new simulated peak. 208 209 while ( 210 abund_min_valley[0] > self.base_line_target 211 or abund_min_valley[1] > self.base_line_target 212 and j <= self.max_interation 213 ): 214 previous_sim_mz, previous_sim_abun = ( 215 previous_peak_obj.gaussian( 216 delta_rp=delta_rp, mz_overlay=self.mz_overlay 217 ) 218 ) 219 220 sim_mz, sim_abun = peak_obj.gaussian( 221 delta_rp=delta_rp, mz_overlay=self.mz_overlay 222 ) 223 224 next_sim_mz, next_sim_abun = next_peak_obj.gaussian( 225 delta_rp=delta_rp, mz_overlay=self.mz_overlay 226 ) 227 228 sim_mz_domain, summed_peaks_abun = self.sum_data( 229 ( 230 (previous_sim_mz, previous_sim_abun), 231 (sim_mz, sim_abun), 232 (next_sim_mz, next_sim_abun), 233 ) 234 ) 235 236 # update_plot(sim_mz_domain, summed_peaks_abun, 0.001) 237 238 # summed_peaks_abun = (sim_abun + next_sim_abun + previous_sim_abun) 239 240 # find appexes location (mz) and magnitude 241 mz_centroid, abund_centroid = self.find_peak_apex( 242 sim_mz_domain, summed_peaks_abun 243 ) 244 245 # find valley location (mz_min_valley) and magnitude (abund_min_valley) 246 summed_peaks_abun = summed_peaks_abun / ( 247 summed_peaks_abun.max() 248 ) 249 mz_min_valley, abund_min_valley = self.find_peak_valley( 250 sim_mz_domain, summed_peaks_abun 251 ) 252 253 if len(abund_min_valley) != 2: 254 break 255 256 delta_rp += self.rp_increments 257 j += 1 258 259 # plot_triplets(mz_centroid,abund_centroid, mz_min_valley, abund_min_valley, sim_mz_domain, summed_peaks_abun ) 260 261 # plot_triplets(mz_centroid,abund_centroid, mz_min_valley, abund_min_valley, sim_mz_domain, summed_peaks_abun ) 262 263 mass_shift_ppp = self.calc_error( 264 mz_centroid[1], peak_obj.mz_exp, 1000000 265 ) 266 # delta_mz = mz_centroid[1] - peak_obj.mz_exp 267 height_shift_per = self.calc_error( 268 abund_centroid[1], peak_obj.abundance, 100 269 ) 270 # excitation_amplitude = str(mass_spectrum_obj.filename.stem).split("ex")[1].split("pc")[0] 271 # ion_time = str(mass_spectrum_obj.filename.stem).split("0pt")[1].split("s")[0] 272 peak_obj.predicted_std = mass_shift_ppp 273 274 results_list.append( 275 { 276 "ms_index_position": peak_obj_idx, 277 "predicted_std": mass_shift_ppp, 278 "mz_exp": peak_obj.mz_exp, 279 "nominal_mz_exp": peak_obj.nominal_mz_exp, 280 "predicted_mz": mz_centroid[1], 281 "s2n": peak_obj.signal_to_noise, 282 "peak_height": peak_obj.abundance, 283 "predicted_peak_height": abund_centroid[1], 284 "peak_height_error": height_shift_per, 285 "resolving_power": peak_obj.resolving_power, 286 # "excitation_amplitude" : excitation_amplitude, 287 # "ion_time" : ion_time 288 } 289 ) 290 291 indexes_without_results.remove(peak_obj_idx) 292 # elif len(mz_centroid) == 3 and len(abund_min_valley) != 2: 293 294 for peak_obj_idx in indexes_without_results: 295 results_list.append( 296 { 297 "ms_index_position": peak_obj_idx, 298 "mz_exp": self.mass_spectrum_obj[peak_obj_idx].mz_exp, 299 "nominal_mz_exp": self.mass_spectrum_obj[ 300 peak_obj_idx 301 ].nominal_mz_exp, 302 "s2n": self.mass_spectrum_obj[peak_obj_idx].signal_to_noise, 303 "peak_height": self.mass_spectrum_obj[peak_obj_idx].abundance, 304 "resolving_power": self.mass_spectrum_obj[ 305 peak_obj_idx 306 ].resolving_power, 307 # "excitation_amplitude" : excitation_amplitude, 308 # "ion_time" : ion_time 309 } 310 ) 311 312 df = DataFrame(results_list).sort_values("mz_exp") 313 314 df.interpolate(method="linear", limit_direction="backward", inplace=True) 315 df.interpolate(method="linear", limit_direction="forward", inplace=True) 316 317 # TODO improve interpolation for missing data 318 # f1 = interpolate.interp1d(x1, y1, kind='quadratic',fill_value="extrapolate") 319 320 for peak_obj_idx in indexes_without_results: 321 predicted_std = df.loc[peak_obj_idx].predicted_std 322 323 self.mass_spectrum_obj[peak_obj_idx].predicted_std = predicted_std 324 325 return df 326 327 def sum_data(self, tuple_mz_abun_list: tuple): 328 """Sum the abundances of the simulated peaks. 329 330 Parameters 331 ------ 332 tuple_mz_abun_list : tuple 333 A tuple containing the mz and abundance lists. 334 335 Returns 336 ------- 337 tuple 338 A tuple containing the summed mz and abundance lists. 339 340 """ 341 all_mz = {} 342 343 for mz_list, abun_list in tuple_mz_abun_list: 344 for index, mz in enumerate(mz_list): 345 abundance = abun_list[index] 346 347 if mz in all_mz: 348 all_mz[mz] = all_mz[mz] + abundance 349 else: 350 all_mz[mz] = abundance 351 352 mz_all = [] 353 abun_all = [] 354 355 for mz in sorted(all_mz): 356 mz_all.append(mz) 357 abun_all.append(all_mz[mz]) 358 359 return array(mz_all), array(abun_all) 360 361 def calc_error(self, mass_ref, mass_sim, factor): 362 """Calculate the error between two values. 363 364 Parameters 365 ---------- 366 mass_ref : float 367 The reference value. 368 mass_sim : float 369 The simulated value. 370 factor : float 371 The factor to multiply the error by. 372 373 Returns 374 ------- 375 float 376 The calculated error. 377 378 379 """ 380 return (mass_sim - mass_ref / mass_ref) * factor 381 382 def find_peak_apex(self, mz, abund): 383 """Find the peak apex. 384 385 Parameters 386 ------ 387 mz : array 388 The mz array. 389 abund : array 390 The abundance array. 391 392 Returns 393 ------- 394 tuple 395 A tuple containing the peak apex mass and abundance. 396 397 """ 398 dy = abund[1:] - abund[:-1] 399 400 # replaces nan for infinity''' 401 indices_nan = where(isnan(abund))[0] 402 403 if indices_nan.size: 404 abund[indices_nan] = inf 405 dy[where(isnan(dy))[0]] = inf 406 407 indexes = where((hstack((dy, 0)) < 0) & (hstack((0, dy)) > 0))[0] 408 409 if indexes.size: 410 return mz[indexes], abund[indexes] 411 412 def find_peak_valley(self, mz, abund): 413 """Find the peak valley. 414 415 Parameters 416 ------ 417 mz : array 418 The mz array. 419 abund : array 420 The abundance array. 421 422 Returns 423 ------- 424 tuple 425 A tuple containing the peak valley mz and abundance. 426 """ 427 dy = abund[1:] - abund[:-1] 428 429 # replaces nan for infinity 430 indices_nan = where(isnan(abund))[0] 431 432 if indices_nan.size: 433 abund[indices_nan] = inf 434 dy[where(isnan(dy))[0]] = inf 435 436 indexes = where((hstack((dy, 0)) > 0) & (hstack((0, dy)) < 0))[0] 437 438 return mz[indexes], abund[indexes]
class
MassErrorPrediction(threading.Thread):
11class MassErrorPrediction(Thread): 12 """Class for mass error prediction. 13 14 Parameters 15 ---------- 16 mass_spectrum : list 17 List of mass spectrum objects. 18 mz_overlay : int, optional 19 The mz overlay value for peak simulation. Default is 10. 20 rp_increments : int, optional 21 The resolving power increments for peak simulation. Default is 10000. 22 base_line_target : float, optional 23 The target value for the baseline resolution. Default is 0.01. 24 max_interation : int, optional 25 The maximum number of iterations for peak simulation. Default is 1000. 26 interpolation : str, optional 27 The interpolation method for missing data. Default is 'linear'. 28 29 Attributes 30 ---------- 31 mass_spectrum_obj : list 32 List of mass spectrum objects. 33 mz_overlay : int 34 The mz overlay value for peak simulation. 35 rp_increments : int 36 The resolving power increments for peak simulation. 37 base_line_target : float 38 The target value for the baseline resolution. 39 max_interation : int 40 The maximum number of iterations for peak simulation. 41 df : DataFrame or None 42 The calculated error distribution dataframe. 43 interpolation : str 44 The interpolation method for missing data. 45 46 Methods 47 ------- 48 * run(). 49 Runs the mass error prediction calculation. 50 * get_results(). 51 Returns the calculated error distribution dataframe. 52 53 """ 54 55 def __init__( 56 self, 57 mass_spectrum, 58 mz_overlay=10, 59 rp_increments=10000, 60 base_line_target: float = 0.01, 61 max_interation=1000, 62 interpolation="linear", 63 ): 64 Thread.__init__(self) 65 66 self.mass_spectrum_obj = mass_spectrum 67 68 self.mz_overlay = mz_overlay 69 70 self.rp_increments = rp_increments 71 72 self.base_line_target = base_line_target 73 74 self.max_interation = max_interation 75 76 self.df = None 77 78 self.interpolation = interpolation 79 80 def run(self): 81 """Runs the mass error prediction calculation.""" 82 self.df = self.calc_error_dist() 83 84 def get_results(self): 85 """Returns the calculated error distribution dataframe.""" 86 87 if not self.df: 88 self.run() 89 90 return self.df 91 92 def calc_error_dist(self): 93 """Calculate the error distribution.""" 94 verbose = self.mass_spectrum_obj.parameters.mass_spectrum.verbose_processing 95 results_list = [] 96 97 indexes_without_results = list(range(len(self.mass_spectrum_obj))) 98 # loop trough mass spectrum 99 for peak_obj_idx, peak_obj in enumerate( 100 tqdm(self.mass_spectrum_obj), disable=not verbose 101 ): 102 # access ms peaks triplets ( peak_obj_idx -1, peak_obj_idx, and peak_obj_idx + 1) 103 # check lower and upper boundaries to not excesses mass spectrum range 104 105 if peak_obj_idx != 0 and peak_obj_idx != len(self.mass_spectrum_obj) - 1: 106 # current peak_obj initialted in the loop expression 107 # geting the peak on the left (previous_peak_obj) and the one in the right position (next_peak_obj) 108 next_peak_obj = self.mass_spectrum_obj[peak_obj_idx + 1] 109 previous_peak_obj = self.mass_spectrum_obj[peak_obj_idx - 1] 110 111 # check mz range defined in max_mz variable and check if peaks have same nominal mz 112 # keeping same mz for better plotting representation only, remove it for production 113 if ( 114 peak_obj.nominal_mz_exp == next_peak_obj.nominal_mz_exp 115 and peak_obj.nominal_mz_exp == previous_peak_obj.nominal_mz_exp 116 ): 117 # simulate peak shape 118 sim_mz, sim_abun = peak_obj.gaussian(mz_overlay=self.mz_overlay) 119 # update_plot(sim_mz,sim_abun, 0.5) 120 121 # simulate peak shape 122 next_sim_mz, next_sim_abun = next_peak_obj.gaussian( 123 mz_overlay=self.mz_overlay 124 ) 125 # update_plot(next_sim_mz, next_sim_abun, 0.5) 126 127 # simulate peak shape 128 previous_sim_mz, previous_sim_abun = previous_peak_obj.gaussian( 129 mz_overlay=self.mz_overlay 130 ) 131 # update_plot(previous_sim_mz, previous_sim_abun, 0.5) 132 133 sim_mz_domain, summed_peaks_abun = self.sum_data( 134 ( 135 (previous_sim_mz, previous_sim_abun), 136 (sim_mz, sim_abun), 137 (next_sim_mz, next_sim_abun), 138 ) 139 ) 140 # update_plot(sim_mz_domain,summed_peaks_abun, 0.5) 141 142 # sum simulated abundances 143 # summed_peaks_abun = (sim_abun + next_sim_abun + previous_sim_abun) 144 145 # normalize abundances to 0-1 146 # summed_peaks_abun = summed_peaks_abun/(max(summed_peaks_abun)) 147 148 # find appexes location (mz) and magnitude 149 mz_centroid, abund_centroid = self.find_peak_apex( 150 sim_mz_domain, summed_peaks_abun 151 ) 152 153 # find valley location (mz_min_valley) and magnitude (abund_min_valley) 154 mz_min_valley, abund_min_valley = self.find_peak_valley( 155 sim_mz_domain, summed_peaks_abun 156 ) 157 158 # clear delta_rp (global implementation) and store choose resolving power increments 159 delta_rp = self.rp_increments 160 161 # used to limited number of iterations 162 i = 0 163 j = 0 164 165 # TODO: fit peak shape and decide best fit #gaussian, lorentz and voigt 166 # plot_triplets(mz_centroid,abund_centroid, mz_min_valley, abund_min_valley, sim_mz_domain, summed_peaks_abun ) 167 if len(mz_centroid) == 2: 168 while len(mz_centroid) < 3 and i <= self.max_interation: 169 previous_sim_mz, previous_sim_abun = ( 170 previous_peak_obj.gaussian( 171 delta_rp=delta_rp, mz_overlay=self.mz_overlay 172 ) 173 ) 174 175 sim_mz, sim_abun = peak_obj.gaussian( 176 delta_rp=delta_rp, mz_overlay=self.mz_overlay 177 ) 178 179 next_sim_mz, next_sim_abun = next_peak_obj.gaussian( 180 delta_rp=delta_rp, mz_overlay=self.mz_overlay 181 ) 182 183 sim_mz_domain, summed_peaks_abun = self.sum_data( 184 ( 185 (previous_sim_mz, previous_sim_abun), 186 (sim_mz, sim_abun), 187 (next_sim_mz, next_sim_abun), 188 ) 189 ) 190 191 # update_plot(sim_mz_domain, summed_peaks_abun, 0.01) 192 193 mz_centroid, abund_centroid = self.find_peak_apex( 194 sim_mz_domain, summed_peaks_abun 195 ) 196 197 delta_rp += self.rp_increments 198 199 i += 1 200 201 mz_min_valley, abund_min_valley = self.find_peak_valley( 202 sim_mz_domain, summed_peaks_abun 203 ) 204 205 if len(mz_centroid) == 3 and len(abund_min_valley) == 2: 206 # increase all three peak resolving power until both valley magnitude is bellow the defined target 207 # calculate peak shapes with the needed resolving power to have a baseline resolution for all peaks 208 # calculate mass difference (ppm) between original centroid and the new simulated peak. 209 210 while ( 211 abund_min_valley[0] > self.base_line_target 212 or abund_min_valley[1] > self.base_line_target 213 and j <= self.max_interation 214 ): 215 previous_sim_mz, previous_sim_abun = ( 216 previous_peak_obj.gaussian( 217 delta_rp=delta_rp, mz_overlay=self.mz_overlay 218 ) 219 ) 220 221 sim_mz, sim_abun = peak_obj.gaussian( 222 delta_rp=delta_rp, mz_overlay=self.mz_overlay 223 ) 224 225 next_sim_mz, next_sim_abun = next_peak_obj.gaussian( 226 delta_rp=delta_rp, mz_overlay=self.mz_overlay 227 ) 228 229 sim_mz_domain, summed_peaks_abun = self.sum_data( 230 ( 231 (previous_sim_mz, previous_sim_abun), 232 (sim_mz, sim_abun), 233 (next_sim_mz, next_sim_abun), 234 ) 235 ) 236 237 # update_plot(sim_mz_domain, summed_peaks_abun, 0.001) 238 239 # summed_peaks_abun = (sim_abun + next_sim_abun + previous_sim_abun) 240 241 # find appexes location (mz) and magnitude 242 mz_centroid, abund_centroid = self.find_peak_apex( 243 sim_mz_domain, summed_peaks_abun 244 ) 245 246 # find valley location (mz_min_valley) and magnitude (abund_min_valley) 247 summed_peaks_abun = summed_peaks_abun / ( 248 summed_peaks_abun.max() 249 ) 250 mz_min_valley, abund_min_valley = self.find_peak_valley( 251 sim_mz_domain, summed_peaks_abun 252 ) 253 254 if len(abund_min_valley) != 2: 255 break 256 257 delta_rp += self.rp_increments 258 j += 1 259 260 # plot_triplets(mz_centroid,abund_centroid, mz_min_valley, abund_min_valley, sim_mz_domain, summed_peaks_abun ) 261 262 # plot_triplets(mz_centroid,abund_centroid, mz_min_valley, abund_min_valley, sim_mz_domain, summed_peaks_abun ) 263 264 mass_shift_ppp = self.calc_error( 265 mz_centroid[1], peak_obj.mz_exp, 1000000 266 ) 267 # delta_mz = mz_centroid[1] - peak_obj.mz_exp 268 height_shift_per = self.calc_error( 269 abund_centroid[1], peak_obj.abundance, 100 270 ) 271 # excitation_amplitude = str(mass_spectrum_obj.filename.stem).split("ex")[1].split("pc")[0] 272 # ion_time = str(mass_spectrum_obj.filename.stem).split("0pt")[1].split("s")[0] 273 peak_obj.predicted_std = mass_shift_ppp 274 275 results_list.append( 276 { 277 "ms_index_position": peak_obj_idx, 278 "predicted_std": mass_shift_ppp, 279 "mz_exp": peak_obj.mz_exp, 280 "nominal_mz_exp": peak_obj.nominal_mz_exp, 281 "predicted_mz": mz_centroid[1], 282 "s2n": peak_obj.signal_to_noise, 283 "peak_height": peak_obj.abundance, 284 "predicted_peak_height": abund_centroid[1], 285 "peak_height_error": height_shift_per, 286 "resolving_power": peak_obj.resolving_power, 287 # "excitation_amplitude" : excitation_amplitude, 288 # "ion_time" : ion_time 289 } 290 ) 291 292 indexes_without_results.remove(peak_obj_idx) 293 # elif len(mz_centroid) == 3 and len(abund_min_valley) != 2: 294 295 for peak_obj_idx in indexes_without_results: 296 results_list.append( 297 { 298 "ms_index_position": peak_obj_idx, 299 "mz_exp": self.mass_spectrum_obj[peak_obj_idx].mz_exp, 300 "nominal_mz_exp": self.mass_spectrum_obj[ 301 peak_obj_idx 302 ].nominal_mz_exp, 303 "s2n": self.mass_spectrum_obj[peak_obj_idx].signal_to_noise, 304 "peak_height": self.mass_spectrum_obj[peak_obj_idx].abundance, 305 "resolving_power": self.mass_spectrum_obj[ 306 peak_obj_idx 307 ].resolving_power, 308 # "excitation_amplitude" : excitation_amplitude, 309 # "ion_time" : ion_time 310 } 311 ) 312 313 df = DataFrame(results_list).sort_values("mz_exp") 314 315 df.interpolate(method="linear", limit_direction="backward", inplace=True) 316 df.interpolate(method="linear", limit_direction="forward", inplace=True) 317 318 # TODO improve interpolation for missing data 319 # f1 = interpolate.interp1d(x1, y1, kind='quadratic',fill_value="extrapolate") 320 321 for peak_obj_idx in indexes_without_results: 322 predicted_std = df.loc[peak_obj_idx].predicted_std 323 324 self.mass_spectrum_obj[peak_obj_idx].predicted_std = predicted_std 325 326 return df 327 328 def sum_data(self, tuple_mz_abun_list: tuple): 329 """Sum the abundances of the simulated peaks. 330 331 Parameters 332 ------ 333 tuple_mz_abun_list : tuple 334 A tuple containing the mz and abundance lists. 335 336 Returns 337 ------- 338 tuple 339 A tuple containing the summed mz and abundance lists. 340 341 """ 342 all_mz = {} 343 344 for mz_list, abun_list in tuple_mz_abun_list: 345 for index, mz in enumerate(mz_list): 346 abundance = abun_list[index] 347 348 if mz in all_mz: 349 all_mz[mz] = all_mz[mz] + abundance 350 else: 351 all_mz[mz] = abundance 352 353 mz_all = [] 354 abun_all = [] 355 356 for mz in sorted(all_mz): 357 mz_all.append(mz) 358 abun_all.append(all_mz[mz]) 359 360 return array(mz_all), array(abun_all) 361 362 def calc_error(self, mass_ref, mass_sim, factor): 363 """Calculate the error between two values. 364 365 Parameters 366 ---------- 367 mass_ref : float 368 The reference value. 369 mass_sim : float 370 The simulated value. 371 factor : float 372 The factor to multiply the error by. 373 374 Returns 375 ------- 376 float 377 The calculated error. 378 379 380 """ 381 return (mass_sim - mass_ref / mass_ref) * factor 382 383 def find_peak_apex(self, mz, abund): 384 """Find the peak apex. 385 386 Parameters 387 ------ 388 mz : array 389 The mz array. 390 abund : array 391 The abundance array. 392 393 Returns 394 ------- 395 tuple 396 A tuple containing the peak apex mass and abundance. 397 398 """ 399 dy = abund[1:] - abund[:-1] 400 401 # replaces nan for infinity''' 402 indices_nan = where(isnan(abund))[0] 403 404 if indices_nan.size: 405 abund[indices_nan] = inf 406 dy[where(isnan(dy))[0]] = inf 407 408 indexes = where((hstack((dy, 0)) < 0) & (hstack((0, dy)) > 0))[0] 409 410 if indexes.size: 411 return mz[indexes], abund[indexes] 412 413 def find_peak_valley(self, mz, abund): 414 """Find the peak valley. 415 416 Parameters 417 ------ 418 mz : array 419 The mz array. 420 abund : array 421 The abundance array. 422 423 Returns 424 ------- 425 tuple 426 A tuple containing the peak valley mz and abundance. 427 """ 428 dy = abund[1:] - abund[:-1] 429 430 # replaces nan for infinity 431 indices_nan = where(isnan(abund))[0] 432 433 if indices_nan.size: 434 abund[indices_nan] = inf 435 dy[where(isnan(dy))[0]] = inf 436 437 indexes = where((hstack((dy, 0)) > 0) & (hstack((0, dy)) < 0))[0] 438 439 return mz[indexes], abund[indexes]
Class for mass error prediction.
Parameters
- mass_spectrum (list): List of mass spectrum objects.
- mz_overlay (int, optional): The mz overlay value for peak simulation. Default is 10.
- rp_increments (int, optional): The resolving power increments for peak simulation. Default is 10000.
- base_line_target (float, optional): The target value for the baseline resolution. Default is 0.01.
- max_interation (int, optional): The maximum number of iterations for peak simulation. Default is 1000.
- interpolation (str, optional): The interpolation method for missing data. Default is 'linear'.
Attributes
- mass_spectrum_obj (list): List of mass spectrum objects.
- mz_overlay (int): The mz overlay value for peak simulation.
- rp_increments (int): The resolving power increments for peak simulation.
- base_line_target (float): The target value for the baseline resolution.
- max_interation (int): The maximum number of iterations for peak simulation.
- df (DataFrame or None): The calculated error distribution dataframe.
- interpolation (str): The interpolation method for missing data.
Methods
- run(). Runs the mass error prediction calculation.
- get_results(). Returns the calculated error distribution dataframe.
MassErrorPrediction( mass_spectrum, mz_overlay=10, rp_increments=10000, base_line_target: float = 0.01, max_interation=1000, interpolation='linear')
55 def __init__( 56 self, 57 mass_spectrum, 58 mz_overlay=10, 59 rp_increments=10000, 60 base_line_target: float = 0.01, 61 max_interation=1000, 62 interpolation="linear", 63 ): 64 Thread.__init__(self) 65 66 self.mass_spectrum_obj = mass_spectrum 67 68 self.mz_overlay = mz_overlay 69 70 self.rp_increments = rp_increments 71 72 self.base_line_target = base_line_target 73 74 self.max_interation = max_interation 75 76 self.df = None 77 78 self.interpolation = interpolation
This constructor should always be called with keyword arguments. Arguments are:
group should be None; reserved for future extension when a ThreadGroup class is implemented.
target is the callable object to be invoked by the run() method. Defaults to None, meaning nothing is called.
name is the thread name. By default, a unique name is constructed of the form "Thread-N" where N is a small decimal number.
args is a list or tuple of arguments for the target invocation. Defaults to ().
kwargs is a dictionary of keyword arguments for the target invocation. Defaults to {}.
If a subclass overrides the constructor, it must make sure to invoke the base class constructor (Thread.__init__()) before doing anything else to the thread.
def
run(self):
80 def run(self): 81 """Runs the mass error prediction calculation.""" 82 self.df = self.calc_error_dist()
Runs the mass error prediction calculation.
def
get_results(self):
84 def get_results(self): 85 """Returns the calculated error distribution dataframe.""" 86 87 if not self.df: 88 self.run() 89 90 return self.df
Returns the calculated error distribution dataframe.
def
calc_error_dist(self):
92 def calc_error_dist(self): 93 """Calculate the error distribution.""" 94 verbose = self.mass_spectrum_obj.parameters.mass_spectrum.verbose_processing 95 results_list = [] 96 97 indexes_without_results = list(range(len(self.mass_spectrum_obj))) 98 # loop trough mass spectrum 99 for peak_obj_idx, peak_obj in enumerate( 100 tqdm(self.mass_spectrum_obj), disable=not verbose 101 ): 102 # access ms peaks triplets ( peak_obj_idx -1, peak_obj_idx, and peak_obj_idx + 1) 103 # check lower and upper boundaries to not excesses mass spectrum range 104 105 if peak_obj_idx != 0 and peak_obj_idx != len(self.mass_spectrum_obj) - 1: 106 # current peak_obj initialted in the loop expression 107 # geting the peak on the left (previous_peak_obj) and the one in the right position (next_peak_obj) 108 next_peak_obj = self.mass_spectrum_obj[peak_obj_idx + 1] 109 previous_peak_obj = self.mass_spectrum_obj[peak_obj_idx - 1] 110 111 # check mz range defined in max_mz variable and check if peaks have same nominal mz 112 # keeping same mz for better plotting representation only, remove it for production 113 if ( 114 peak_obj.nominal_mz_exp == next_peak_obj.nominal_mz_exp 115 and peak_obj.nominal_mz_exp == previous_peak_obj.nominal_mz_exp 116 ): 117 # simulate peak shape 118 sim_mz, sim_abun = peak_obj.gaussian(mz_overlay=self.mz_overlay) 119 # update_plot(sim_mz,sim_abun, 0.5) 120 121 # simulate peak shape 122 next_sim_mz, next_sim_abun = next_peak_obj.gaussian( 123 mz_overlay=self.mz_overlay 124 ) 125 # update_plot(next_sim_mz, next_sim_abun, 0.5) 126 127 # simulate peak shape 128 previous_sim_mz, previous_sim_abun = previous_peak_obj.gaussian( 129 mz_overlay=self.mz_overlay 130 ) 131 # update_plot(previous_sim_mz, previous_sim_abun, 0.5) 132 133 sim_mz_domain, summed_peaks_abun = self.sum_data( 134 ( 135 (previous_sim_mz, previous_sim_abun), 136 (sim_mz, sim_abun), 137 (next_sim_mz, next_sim_abun), 138 ) 139 ) 140 # update_plot(sim_mz_domain,summed_peaks_abun, 0.5) 141 142 # sum simulated abundances 143 # summed_peaks_abun = (sim_abun + next_sim_abun + previous_sim_abun) 144 145 # normalize abundances to 0-1 146 # summed_peaks_abun = summed_peaks_abun/(max(summed_peaks_abun)) 147 148 # find appexes location (mz) and magnitude 149 mz_centroid, abund_centroid = self.find_peak_apex( 150 sim_mz_domain, summed_peaks_abun 151 ) 152 153 # find valley location (mz_min_valley) and magnitude (abund_min_valley) 154 mz_min_valley, abund_min_valley = self.find_peak_valley( 155 sim_mz_domain, summed_peaks_abun 156 ) 157 158 # clear delta_rp (global implementation) and store choose resolving power increments 159 delta_rp = self.rp_increments 160 161 # used to limited number of iterations 162 i = 0 163 j = 0 164 165 # TODO: fit peak shape and decide best fit #gaussian, lorentz and voigt 166 # plot_triplets(mz_centroid,abund_centroid, mz_min_valley, abund_min_valley, sim_mz_domain, summed_peaks_abun ) 167 if len(mz_centroid) == 2: 168 while len(mz_centroid) < 3 and i <= self.max_interation: 169 previous_sim_mz, previous_sim_abun = ( 170 previous_peak_obj.gaussian( 171 delta_rp=delta_rp, mz_overlay=self.mz_overlay 172 ) 173 ) 174 175 sim_mz, sim_abun = peak_obj.gaussian( 176 delta_rp=delta_rp, mz_overlay=self.mz_overlay 177 ) 178 179 next_sim_mz, next_sim_abun = next_peak_obj.gaussian( 180 delta_rp=delta_rp, mz_overlay=self.mz_overlay 181 ) 182 183 sim_mz_domain, summed_peaks_abun = self.sum_data( 184 ( 185 (previous_sim_mz, previous_sim_abun), 186 (sim_mz, sim_abun), 187 (next_sim_mz, next_sim_abun), 188 ) 189 ) 190 191 # update_plot(sim_mz_domain, summed_peaks_abun, 0.01) 192 193 mz_centroid, abund_centroid = self.find_peak_apex( 194 sim_mz_domain, summed_peaks_abun 195 ) 196 197 delta_rp += self.rp_increments 198 199 i += 1 200 201 mz_min_valley, abund_min_valley = self.find_peak_valley( 202 sim_mz_domain, summed_peaks_abun 203 ) 204 205 if len(mz_centroid) == 3 and len(abund_min_valley) == 2: 206 # increase all three peak resolving power until both valley magnitude is bellow the defined target 207 # calculate peak shapes with the needed resolving power to have a baseline resolution for all peaks 208 # calculate mass difference (ppm) between original centroid and the new simulated peak. 209 210 while ( 211 abund_min_valley[0] > self.base_line_target 212 or abund_min_valley[1] > self.base_line_target 213 and j <= self.max_interation 214 ): 215 previous_sim_mz, previous_sim_abun = ( 216 previous_peak_obj.gaussian( 217 delta_rp=delta_rp, mz_overlay=self.mz_overlay 218 ) 219 ) 220 221 sim_mz, sim_abun = peak_obj.gaussian( 222 delta_rp=delta_rp, mz_overlay=self.mz_overlay 223 ) 224 225 next_sim_mz, next_sim_abun = next_peak_obj.gaussian( 226 delta_rp=delta_rp, mz_overlay=self.mz_overlay 227 ) 228 229 sim_mz_domain, summed_peaks_abun = self.sum_data( 230 ( 231 (previous_sim_mz, previous_sim_abun), 232 (sim_mz, sim_abun), 233 (next_sim_mz, next_sim_abun), 234 ) 235 ) 236 237 # update_plot(sim_mz_domain, summed_peaks_abun, 0.001) 238 239 # summed_peaks_abun = (sim_abun + next_sim_abun + previous_sim_abun) 240 241 # find appexes location (mz) and magnitude 242 mz_centroid, abund_centroid = self.find_peak_apex( 243 sim_mz_domain, summed_peaks_abun 244 ) 245 246 # find valley location (mz_min_valley) and magnitude (abund_min_valley) 247 summed_peaks_abun = summed_peaks_abun / ( 248 summed_peaks_abun.max() 249 ) 250 mz_min_valley, abund_min_valley = self.find_peak_valley( 251 sim_mz_domain, summed_peaks_abun 252 ) 253 254 if len(abund_min_valley) != 2: 255 break 256 257 delta_rp += self.rp_increments 258 j += 1 259 260 # plot_triplets(mz_centroid,abund_centroid, mz_min_valley, abund_min_valley, sim_mz_domain, summed_peaks_abun ) 261 262 # plot_triplets(mz_centroid,abund_centroid, mz_min_valley, abund_min_valley, sim_mz_domain, summed_peaks_abun ) 263 264 mass_shift_ppp = self.calc_error( 265 mz_centroid[1], peak_obj.mz_exp, 1000000 266 ) 267 # delta_mz = mz_centroid[1] - peak_obj.mz_exp 268 height_shift_per = self.calc_error( 269 abund_centroid[1], peak_obj.abundance, 100 270 ) 271 # excitation_amplitude = str(mass_spectrum_obj.filename.stem).split("ex")[1].split("pc")[0] 272 # ion_time = str(mass_spectrum_obj.filename.stem).split("0pt")[1].split("s")[0] 273 peak_obj.predicted_std = mass_shift_ppp 274 275 results_list.append( 276 { 277 "ms_index_position": peak_obj_idx, 278 "predicted_std": mass_shift_ppp, 279 "mz_exp": peak_obj.mz_exp, 280 "nominal_mz_exp": peak_obj.nominal_mz_exp, 281 "predicted_mz": mz_centroid[1], 282 "s2n": peak_obj.signal_to_noise, 283 "peak_height": peak_obj.abundance, 284 "predicted_peak_height": abund_centroid[1], 285 "peak_height_error": height_shift_per, 286 "resolving_power": peak_obj.resolving_power, 287 # "excitation_amplitude" : excitation_amplitude, 288 # "ion_time" : ion_time 289 } 290 ) 291 292 indexes_without_results.remove(peak_obj_idx) 293 # elif len(mz_centroid) == 3 and len(abund_min_valley) != 2: 294 295 for peak_obj_idx in indexes_without_results: 296 results_list.append( 297 { 298 "ms_index_position": peak_obj_idx, 299 "mz_exp": self.mass_spectrum_obj[peak_obj_idx].mz_exp, 300 "nominal_mz_exp": self.mass_spectrum_obj[ 301 peak_obj_idx 302 ].nominal_mz_exp, 303 "s2n": self.mass_spectrum_obj[peak_obj_idx].signal_to_noise, 304 "peak_height": self.mass_spectrum_obj[peak_obj_idx].abundance, 305 "resolving_power": self.mass_spectrum_obj[ 306 peak_obj_idx 307 ].resolving_power, 308 # "excitation_amplitude" : excitation_amplitude, 309 # "ion_time" : ion_time 310 } 311 ) 312 313 df = DataFrame(results_list).sort_values("mz_exp") 314 315 df.interpolate(method="linear", limit_direction="backward", inplace=True) 316 df.interpolate(method="linear", limit_direction="forward", inplace=True) 317 318 # TODO improve interpolation for missing data 319 # f1 = interpolate.interp1d(x1, y1, kind='quadratic',fill_value="extrapolate") 320 321 for peak_obj_idx in indexes_without_results: 322 predicted_std = df.loc[peak_obj_idx].predicted_std 323 324 self.mass_spectrum_obj[peak_obj_idx].predicted_std = predicted_std 325 326 return df
Calculate the error distribution.
def
sum_data(self, tuple_mz_abun_list: tuple):
328 def sum_data(self, tuple_mz_abun_list: tuple): 329 """Sum the abundances of the simulated peaks. 330 331 Parameters 332 ------ 333 tuple_mz_abun_list : tuple 334 A tuple containing the mz and abundance lists. 335 336 Returns 337 ------- 338 tuple 339 A tuple containing the summed mz and abundance lists. 340 341 """ 342 all_mz = {} 343 344 for mz_list, abun_list in tuple_mz_abun_list: 345 for index, mz in enumerate(mz_list): 346 abundance = abun_list[index] 347 348 if mz in all_mz: 349 all_mz[mz] = all_mz[mz] + abundance 350 else: 351 all_mz[mz] = abundance 352 353 mz_all = [] 354 abun_all = [] 355 356 for mz in sorted(all_mz): 357 mz_all.append(mz) 358 abun_all.append(all_mz[mz]) 359 360 return array(mz_all), array(abun_all)
Sum the abundances of the simulated peaks.
Parameters
- tuple_mz_abun_list (tuple): A tuple containing the mz and abundance lists.
Returns
- tuple: A tuple containing the summed mz and abundance lists.
def
calc_error(self, mass_ref, mass_sim, factor):
362 def calc_error(self, mass_ref, mass_sim, factor): 363 """Calculate the error between two values. 364 365 Parameters 366 ---------- 367 mass_ref : float 368 The reference value. 369 mass_sim : float 370 The simulated value. 371 factor : float 372 The factor to multiply the error by. 373 374 Returns 375 ------- 376 float 377 The calculated error. 378 379 380 """ 381 return (mass_sim - mass_ref / mass_ref) * factor
Calculate the error between two values.
Parameters
- mass_ref (float): The reference value.
- mass_sim (float): The simulated value.
- factor (float): The factor to multiply the error by.
Returns
- float: The calculated error.
def
find_peak_apex(self, mz, abund):
383 def find_peak_apex(self, mz, abund): 384 """Find the peak apex. 385 386 Parameters 387 ------ 388 mz : array 389 The mz array. 390 abund : array 391 The abundance array. 392 393 Returns 394 ------- 395 tuple 396 A tuple containing the peak apex mass and abundance. 397 398 """ 399 dy = abund[1:] - abund[:-1] 400 401 # replaces nan for infinity''' 402 indices_nan = where(isnan(abund))[0] 403 404 if indices_nan.size: 405 abund[indices_nan] = inf 406 dy[where(isnan(dy))[0]] = inf 407 408 indexes = where((hstack((dy, 0)) < 0) & (hstack((0, dy)) > 0))[0] 409 410 if indexes.size: 411 return mz[indexes], abund[indexes]
Find the peak apex.
Parameters
- mz (array): The mz array.
- abund (array): The abundance array.
Returns
- tuple: A tuple containing the peak apex mass and abundance.
def
find_peak_valley(self, mz, abund):
413 def find_peak_valley(self, mz, abund): 414 """Find the peak valley. 415 416 Parameters 417 ------ 418 mz : array 419 The mz array. 420 abund : array 421 The abundance array. 422 423 Returns 424 ------- 425 tuple 426 A tuple containing the peak valley mz and abundance. 427 """ 428 dy = abund[1:] - abund[:-1] 429 430 # replaces nan for infinity 431 indices_nan = where(isnan(abund))[0] 432 433 if indices_nan.size: 434 abund[indices_nan] = inf 435 dy[where(isnan(dy))[0]] = inf 436 437 indexes = where((hstack((dy, 0)) > 0) & (hstack((0, dy)) < 0))[0] 438 439 return mz[indexes], abund[indexes]
Find the peak valley.
Parameters
- mz (array): The mz array.
- abund (array): The abundance array.
Returns
- tuple: A tuple containing the peak valley mz and abundance.