Package specific_fns :: Package model_free :: Module bmrb
[hide private]
[frames] | no frames]

Source Code for Module specific_fns.model_free.bmrb

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2009-2013 Edward d'Auvergne                                   # 
  4  #                                                                             # 
  5  # This file is part of the program relax (http://www.nmr-relax.com).          # 
  6  #                                                                             # 
  7  # This program is free software: you can redistribute it and/or modify        # 
  8  # it under the terms of the GNU General Public License as published by        # 
  9  # the Free Software Foundation, either version 3 of the License, or           # 
 10  # (at your option) any later version.                                         # 
 11  #                                                                             # 
 12  # This program is distributed in the hope that it will be useful,             # 
 13  # but WITHOUT ANY WARRANTY; without even the implied warranty of              # 
 14  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
 15  # GNU General Public License for more details.                                # 
 16  #                                                                             # 
 17  # You should have received a copy of the GNU General Public License           # 
 18  # along with this program.  If not, see <http://www.gnu.org/licenses/>.       # 
 19  #                                                                             # 
 20  ############################################################################### 
 21   
 22  # Module docstring. 
 23  """The BMRB methods of the specific API for model-free analysis.""" 
 24   
 25  # Dependency check module. 
 26  import dep_check 
 27   
 28  # Python module imports. 
 29  from math import pi 
 30  from numpy import int32, zeros 
 31  import string 
 32  from warnings import warn 
 33   
 34  # relax module imports. 
 35  if dep_check.bmrblib_module: 
 36      import bmrblib 
 37  from generic_fns import bmrb, diffusion_tensor, exp_info, mol_res_spin, pipes, relax_data 
 38  from generic_fns.interatomic import return_interatom_list 
 39  from generic_fns.mol_res_spin import get_molecule_names, spin_loop 
 40  from relax_errors import RelaxError 
 41  from relax_warnings import RelaxWarning 
 42   
 43   
44 -class Bmrb:
45 """Class containing methods related to BMRB STAR file reading and writing.""" 46
47 - def _from_bmrb_model(self, name=None):
48 """The model-free model name to BMRB name mapping. 49 50 @keyword name: The BMRB model name. 51 @type name: str 52 @return: The corresponding model-free model name. 53 @rtype: str 54 """ 55 56 # The same name. 57 if name in ['m0', 'm1', 'm2', 'm3', 'm4', 'm5', 'm6', 'm7', 'm8', 'm9', 'm0', 'tm1', 'tm2', 'tm3', 'tm4', 'tm5', 'tm6', 'tm7', 'tm8', 'tm9']: 58 return name 59 60 # Conversion of Modelfree4 (and relax) model numbers. 61 if name in ['0', '1', '2', '3', '4', '5', '6', '7', '8', '9']: 62 return 'm' + name 63 64 # The BMRB to model-free model name map. 65 map = {'': 'm0', 66 'S2': 'm1', 67 'S2, te': 'm2', 68 'S2, Rex': 'm3', 69 'S2, te, Rex': 'm4', 70 'S2, te, S2f': 'm5', 71 'S2f, S2, ts': 'm5', 72 'S2f, tf, S2, ts': 'm6', 73 'S2f, S2, ts, Rex': 'm7', 74 'S2f, tf, S2, ts, Rex': 'm8', 75 'Rex': 'm9', 76 'tm': 'tm0', 77 'tm, S2': 'tm1', 78 'tm, S2, te': 'tm2', 79 'tm, S2, Rex': 'tm3', 80 'tm, S2, te, Rex': 'tm4', 81 'tm, S2f, S2, ts': 'tm5', 82 'tm, S2f, tf, S2, ts': 'tm6', 83 'tm, S2f, S2, ts, Rex': 'tm7', 84 'tm, S2f, tf, S2, ts, Rex': 'tm8', 85 'tm, Rex': 'tm9' 86 } 87 88 # Loop over the dictionary. 89 for item in map.items(): 90 # Normal match. 91 if item[0] == name: 92 return item[1] 93 94 # No whitespace. 95 if item[0].replace(' ', '') == name: 96 return item[1] 97 98 # Should not be here! 99 if name: 100 raise RelaxError("The BMRB model-free model '%s' is unknown." % name)
101 102
103 - def _sf_model_free_read(self, star, sample_conditions=None):
104 """Fill the spin containers with the model-free data from the saveframe records. 105 106 @param star: The NMR-STAR dictionary object. 107 @type star: NMR_STAR instance 108 @keyword sample_conditions: The sample condition label to read. Only one sample condition can be read per data pipe. 109 @type sample_conditions: None or str 110 """ 111 112 # The list of model-free parameters (both bmrblib names and relax names). 113 mf_bmrb_key = ['bond_length', 'local_tm', 's2', 's2f', 's2s', 'te', 'tf', 'ts', 'rex', 'chi2'] 114 mf_params = ['r', 'local_tm', 's2', 's2f', 's2s', 'te', 'tf', 'ts', 'rex', 'chi2'] 115 116 # Get the entities. 117 for data in star.model_free.loop(): 118 # Store the keys. 119 keys = data.keys() 120 121 # Sample conditions do not match (remove the $ sign). 122 if 'sample_cond_list_label' in keys and sample_conditions and data['sample_cond_list_label'].replace('$', '') != sample_conditions: 123 continue 124 125 # Global data. 126 if 'global_chi2' in keys: 127 setattr(cdp, 'chi2', data['global_chi2']) 128 129 # The number of spins. 130 N = bmrb.num_spins(data) 131 132 # No data in the saveframe. 133 if N == 0: 134 continue 135 136 # The molecule names. 137 mol_names = bmrb.molecule_names(data, N) 138 139 # Missing atom names. 140 if 'atom_names' not in keys or data['atom_names'] == None: 141 data['atom_names'] = [None] * N 142 143 # Generate the sequence if needed. 144 bmrb.generate_sequence(N, spin_names=data['atom_names'], res_nums=data['res_nums'], res_names=data['res_names'], mol_names=mol_names) 145 146 # Correlation time scaling. 147 table = {'s': 1.0, 148 'ns': 1e-9, 149 'ps': 1e-12} 150 te_scale = 1.0 151 if data['te_units']: 152 te_scale = table[data['te_units']] 153 154 # Fast correlation time scaling. 155 if data['tf_units']: 156 tf_scale = table[data['tf_units']] 157 else: 158 tf_scale = te_scale 159 160 # Slow correlation time scaling. 161 if data['ts_units']: 162 ts_scale = table[data['ts_units']] 163 else: 164 ts_scale = te_scale 165 166 # Rex scaling. 167 rex_scale = 1.0 168 if hasattr(cdp, 'frq') and len(cdp.frq): 169 rex_scale = 1.0 / (2.0*pi*cdp.frq[cdp.ri_ids[0]])**2 170 171 # Loop over the spins. 172 for i in range(N): 173 # Generate a spin ID. 174 spin_id = mol_res_spin.generate_spin_id_unique(mol_name=mol_names[i], res_name=data['res_names'][i], res_num=data['res_nums'][i], spin_name=data['atom_names'][i]) 175 176 # Obtain the spin. 177 spin = mol_res_spin.return_spin(spin_id) 178 179 # No spin?!? 180 if spin == None: 181 raise RelaxError("The spin '%s' does not exist." % spin_id) 182 183 # Loop over and set the model-free parameters. 184 for j in range(len(mf_params)): 185 # The parameter. 186 param = mf_params[j] 187 188 # No parameter. 189 if not mf_bmrb_key[j] in keys: 190 continue 191 192 # The parameter and its value. 193 if data[mf_bmrb_key[j]] != None: 194 # The value. 195 value = data[mf_bmrb_key[j]][i] 196 197 # Parameter scaling. 198 if value != None: 199 if param == 'te': 200 value = value * te_scale 201 elif param == 'tf': 202 value = value * tf_scale 203 elif param == 'ts': 204 value = value * ts_scale 205 elif param == 'rex': 206 value = value * rex_scale 207 208 # No parameter value. 209 else: 210 value = None 211 212 # Set the parameter. 213 setattr(spin, param, value) 214 215 # The error. 216 mf_bmrb_key_err = mf_bmrb_key[j] + '_err' 217 error = None 218 if mf_bmrb_key_err in keys and data[mf_bmrb_key_err] != None: 219 error = data[mf_bmrb_key_err][i] 220 221 # Error scaling. 222 if error != None: 223 if param == 'te': 224 error = error * te_scale 225 elif param == 'tf': 226 error = error * tf_scale 227 elif param == 'ts': 228 error = error * ts_scale 229 elif param == 'rex': 230 error = error * rex_scale 231 232 # Set the error. 233 mf_param_err = param + '_err' 234 if mf_bmrb_key_err in keys and data[mf_bmrb_key_err] != None: 235 setattr(spin, mf_param_err, error) 236 237 # The model. 238 if data['model_fit'] != None and data['model_fit'][i] != None: 239 model = self._from_bmrb_model(data['model_fit'][i]) 240 setattr(spin, 'model', model) 241 242 # The equation and parameters. 243 equation, params = self._model_map(model) 244 setattr(spin, 'equation', equation) 245 setattr(spin, 'params', params) 246 247 # Convert te values which should be ts! 248 if hasattr(spin, 'model') and spin.model in ['m5', 'm6', 'm7', 'm8'] and hasattr(spin, 'te') and spin.te != None: 249 # Change the parameter name of te to ts. 250 spin.ts = spin.te 251 if hasattr(spin, 'te_err'): 252 spin.ts_err = spin.te_err 253 254 # Set the te and te_err values to None. 255 spin.te = None 256 spin.te_err = None 257 258 # The element. 259 if'atom_types' in keys and data['atom_types'] != None: 260 setattr(spin, 'element', data['atom_types'][i]) 261 262 # Heteronucleus type. 263 if'atom_types' in keys and data['atom_types'] != None and data['atom_types'][i] != None and 'isotope' in keys and data['isotope'] != None: 264 # The isotope number. 265 iso_num = data['isotope'][i] 266 267 # No isotope number. 268 iso_table = {'C': 13, 'N': 15} 269 if not data['isotope'][i]: 270 iso_num = iso_table[data['atom_types'][i]] 271 272 # Set the type. 273 setattr(spin, 'isotope', str(iso_num) + data['atom_types'][i])
274 275
276 - def _sf_csa_read(self, star):
277 """Place the CSA data from the saveframe records into the spin container. 278 279 @param star: The NMR-STAR dictionary object. 280 @type star: NMR_STAR instance 281 """ 282 283 # Get the entities. 284 for data in star.chem_shift_anisotropy.loop(): 285 # The number of spins. 286 N = bmrb.num_spins(data) 287 288 # No data in the saveframe. 289 if N == 0: 290 continue 291 292 # The molecule names. 293 mol_names = bmrb.molecule_names(data, N) 294 295 # Loop over the spins. 296 for i in range(len(data['data_ids'])): 297 # Generate a spin ID. 298 spin_id = mol_res_spin.generate_spin_id_unique(mol_name=mol_names[i], res_num=data['res_nums'][i], spin_name=data['atom_names'][i]) 299 300 # Obtain the spin. 301 spin = mol_res_spin.return_spin(spin_id) 302 303 # The CSA value (converted from ppm). 304 setattr(spin, 'csa', data['csa'][i] * 1e-6)
305 306
307 - def _to_bmrb_model(self, name=None):
308 """Convert the model-free model name to the BMRB name. 309 310 @keyword name: The model-free model name. 311 @type name: str 312 @return: The corresponding BMRB model name. 313 @rtype: str 314 """ 315 316 # The relax to BMRB model-free model name map. 317 map = {'m0': '', 318 'm1': 'S2', 319 'm2': 'S2, te', 320 'm3': 'S2, Rex', 321 'm4': 'S2, te, Rex', 322 'm5': 'S2f, S2, ts', 323 'm6': 'S2f, tf, S2, ts', 324 'm7': 'S2f, S2, ts, Rex', 325 'm8': 'S2f, tf, S2, ts, Rex', 326 'm9': 'Rex', 327 'tm0': 'tm', 328 'tm1': 'tm, S2', 329 'tm2': 'tm, S2, te', 330 'tm3': 'tm, S2, Rex', 331 'tm4': 'tm, S2, te, Rex', 332 'tm5': 'tm, S2f, S2, ts', 333 'tm6': 'tm, S2f, tf, S2, ts', 334 'tm7': 'tm, S2f, S2, ts, Rex', 335 'tm8': 'tm, S2f, tf, S2, ts, Rex', 336 'tm9': 'tm, Rex' 337 } 338 339 # No match. 340 if name not in map.keys(): 341 raise RelaxError("The model-free model '%s' is unknown." % name) 342 343 # Return the BMRB model name. 344 return map[name]
345 346
347 - def bmrb_read(self, file_path, version=None, sample_conditions=None):
348 """Read the model-free results from a BMRB NMR-STAR v3.1 formatted file. 349 350 @param file_path: The full file path. 351 @type file_path: str 352 @keyword version: The BMRB version to force the reading. 353 @type version: None or str 354 @keyword sample_conditions: The sample condition label to read. Only one sample condition can be read per data pipe. 355 @type sample_conditions: None or str 356 """ 357 358 # Initialise the NMR-STAR data object. 359 star = bmrblib.create_nmr_star('relax_model_free_results', file_path, version) 360 361 # Read the contents of the STAR formatted file. 362 star.read() 363 364 # The sample conditions. 365 sample_conds = bmrb.list_sample_conditions(star) 366 if sample_conditions and sample_conditions not in sample_conds: 367 raise RelaxError("The sample conditions label '%s' does not correspond to any of the labels in the file: %s" % (sample_conditions, sample_conds)) 368 if not sample_conditions and len(sample_conds) > 1: 369 raise RelaxError("Only one of the sample conditions in %s can be loaded per relax data pipe." % sample_conds) 370 371 # The diffusion tensor. 372 diffusion_tensor.bmrb_read(star) 373 374 # Generate the molecule and residue containers from the entity records. 375 mol_res_spin.bmrb_read(star) 376 377 # Read the relaxation data saveframes. 378 relax_data.bmrb_read(star, sample_conditions=sample_conditions) 379 380 # Read the model-free data saveframes. 381 self._sf_model_free_read(star, sample_conditions=sample_conditions) 382 383 # Read the CSA data saveframes. 384 self._sf_csa_read(star)
385 386
387 - def bmrb_write(self, file_path, version=None):
388 """Write the model-free results to a BMRB NMR-STAR v3.1 formatted file. 389 390 @param file_path: The full file path. 391 @type file_path: str 392 @keyword version: The BMRB NMR-STAR dictionary format to output to. 393 @type version: str 394 """ 395 396 # Alias the current data pipe. 397 cdp = pipes.get_pipe() 398 399 # Initialise the NMR-STAR data object. 400 star = bmrblib.create_nmr_star('relax_model_free_results', file_path, version) 401 402 # Global minimisation stats. 403 global_chi2 = None 404 if hasattr(cdp, 'chi2'): 405 global_chi2 = cdp.chi2 406 407 # Rex frq. 408 rex_frq = cdp.frq[cdp.ri_ids[0]] 409 410 # Initialise the spin specific data lists. 411 mol_name_list = [] 412 res_num_list = [] 413 res_name_list = [] 414 atom_name_list = [] 415 416 csa_list = [] 417 r_list = [] 418 isotope_list = [] 419 element_list = [] 420 421 local_tm_list = [] 422 s2_list = [] 423 s2f_list = [] 424 s2s_list = [] 425 te_list = [] 426 tf_list = [] 427 ts_list = [] 428 rex_list = [] 429 430 local_tm_err_list = [] 431 s2_err_list = [] 432 s2f_err_list = [] 433 s2s_err_list = [] 434 te_err_list = [] 435 tf_err_list = [] 436 ts_err_list = [] 437 rex_err_list = [] 438 439 chi2_list = [] 440 model_list = [] 441 442 # Store the spin specific data in lists for later use. 443 for spin, mol_name, res_num, res_name, spin_id in spin_loop(full_info=True, return_id=True): 444 # Skip the protons. 445 if spin.name == 'H' or (hasattr(spin, 'element') and spin.element == 'H'): 446 warn(RelaxWarning("Skipping the proton spin '%s'." % spin_id)) 447 continue 448 449 # Check the data for None (not allowed in BMRB!). 450 if res_num == None: 451 raise RelaxError("For the BMRB, the residue of spin '%s' must be numbered." % spin_id) 452 if res_name == None: 453 raise RelaxError("For the BMRB, the residue of spin '%s' must be named." % spin_id) 454 if spin.name == None: 455 raise RelaxError("For the BMRB, the spin '%s' must be named." % spin_id) 456 if not hasattr(spin, 'isotope') or spin.isotope == None: 457 raise RelaxError("For the BMRB, the spin isotope type of '%s' must be specified." % spin_id) 458 if not hasattr(spin, 'element') or spin.element == None: 459 raise RelaxError("For the BMRB, the spin element type of '%s' must be specified. Please use the spin user function for setting the element type." % spin_id) 460 461 # The molecule/residue/spin info. 462 mol_name_list.append(mol_name) 463 res_num_list.append(res_num) 464 res_name_list.append(res_name) 465 atom_name_list.append(spin.name) 466 467 # CSA values. 468 if hasattr(spin, 'csa'): 469 csa_list.append(spin.csa * 1e6) # In ppm. 470 else: 471 csa_list.append(None) 472 473 # Interatomic distances. 474 interatoms = return_interatom_list(spin_id) 475 for i in range(len(interatoms)): 476 # No relaxation mechanism. 477 if not interatoms[i].dipole_pair: 478 continue 479 480 # Add the interatomic distance. 481 if hasattr(interatoms[i], 'r'): 482 r_list.append(interatoms[i].r) 483 else: 484 r_list.append(None) 485 486 # Stop adding. 487 break 488 489 # The nuclear isotope. 490 if hasattr(spin, 'isotope'): 491 isotope_list.append(int(spin.isotope.strip(string.ascii_letters))) 492 else: 493 isotope_list.append(None) 494 495 # The element. 496 if hasattr(spin, 'element'): 497 element_list.append(spin.element) 498 else: 499 element_list.append(None) 500 501 # Diffusion tensor. 502 if hasattr(spin, 'local_tm'): 503 local_tm_list.append(spin.local_tm) 504 else: 505 local_tm_list.append(None) 506 if hasattr(spin, 'local_tm_err'): 507 local_tm_err_list.append(spin.local_tm_err) 508 else: 509 local_tm_err_list.append(None) 510 511 # Model-free parameter values. 512 s2_list.append(spin.s2) 513 s2f_list.append(spin.s2f) 514 s2s_list.append(spin.s2s) 515 te_list.append(spin.te) 516 tf_list.append(spin.tf) 517 ts_list.append(spin.ts) 518 if spin.rex == None: 519 rex_list.append(None) 520 else: 521 rex_list.append(spin.rex / (2.0*pi*rex_frq**2)) 522 523 # Model-free parameter errors. 524 params = ['s2', 's2f', 's2s', 'te', 'tf', 'ts', 'rex'] 525 for param in params: 526 # The error list. 527 err_list = locals()[param+'_err_list'] 528 529 # Append the error. 530 if hasattr(spin, param+'_err'): 531 # The value. 532 val = getattr(spin, param+'_err') 533 534 # Scaling. 535 if param == 'rex' and val != None: 536 val = val / (2.0*pi*rex_frq**2) 537 538 # Append. 539 err_list.append(val) 540 541 # Or None. 542 else: 543 err_list.append(None) 544 545 546 # Opt stats. 547 chi2_list.append(spin.chi2) 548 549 # Model-free model. 550 model_list.append(self._to_bmrb_model(spin.model)) 551 552 # Convert the molecule names into the entity IDs. 553 entity_ids = zeros(len(mol_name_list), int32) 554 mol_names = get_molecule_names() 555 for i in range(len(mol_name_list)): 556 for j in range(len(mol_names)): 557 if mol_name_list[i] == mol_names[j]: 558 entity_ids[i] = j+1 559 560 561 # Create Supergroup 2 : The citations. 562 ###################################### 563 564 # Generate the citations saveframe. 565 exp_info.bmrb_write_citations(star) 566 567 568 # Create Supergroup 3 : The molecular assembly saveframes. 569 ########################################################## 570 571 # Generate the entity saveframe. 572 mol_res_spin.bmrb_write_entity(star) 573 574 575 # Create Supergroup 4: The experimental descriptions saveframes. 576 ################################################################# 577 578 # Generate the method saveframes. 579 exp_info.bmrb_write_methods(star) 580 581 # Generate the software saveframe. 582 software_ids, software_labels = exp_info.bmrb_write_software(star) 583 584 585 # Create Supergroup 5 : The NMR parameters saveframes. 586 ###################################################### 587 588 # Generate the CSA saveframe. 589 star.chem_shift_anisotropy.add(entity_ids=entity_ids, res_nums=res_num_list, res_names=res_name_list, atom_names=atom_name_list, atom_types=element_list, isotope=isotope_list, csa=csa_list) 590 591 592 # Create Supergroup 6 : The kinetic data saveframes. 593 #################################################### 594 595 # Generate the relaxation data saveframes. 596 relax_data.bmrb_write(star) 597 598 599 # Create Supergroup 7 : The thermodynamics saveframes. 600 ###################################################### 601 602 # Get the relax software id. 603 for i in range(len(software_ids)): 604 if software_labels[i] == 'relax': 605 software_id = software_ids[i] 606 607 # Generate the model-free data saveframe. 608 star.model_free.add(global_chi2=global_chi2, software_ids=[software_id], software_labels=['relax'], entity_ids=entity_ids, res_nums=res_num_list, res_names=res_name_list, atom_names=atom_name_list, atom_types=element_list, isotope=isotope_list, bond_length=r_list, local_tc=local_tm_list, s2=s2_list, s2f=s2f_list, s2s=s2s_list, te=te_list, tf=tf_list, ts=ts_list, rex=rex_list, local_tc_err=local_tm_err_list, s2_err=s2_err_list, s2f_err=s2f_err_list, s2s_err=s2s_err_list, te_err=te_err_list, tf_err=tf_err_list, ts_err=ts_err_list, rex_err=rex_err_list, rex_frq=rex_frq, chi2=chi2_list, model_fit=model_list) 609 610 611 # Create Supergroup 8 : The structure determination saveframes. 612 ############################################################### 613 614 # Generate the diffusion tensor saveframes. 615 diffusion_tensor.bmrb_write(star) 616 617 618 # Write the contents to the STAR formatted file. 619 star.write()
620