pipe_control.dasha

1 ############################################################################### 2 # # 3 # Copyright (C) 2001-2009,2011-2012 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 """Module for interfacing with Dasha.""" 24 25 # Dependencies. 26 import dep_check 27 28 # Python module imports. 29 from math import pi 30 from os import F_OK, access, chdir, getcwd, sep 31 PIPE, Popen = None, None 32 if dep_check.subprocess_module: 33 from subprocess import PIPE, Popen 34 import sys 35 36 # relax module imports. 37 from lib.errors import RelaxDirError, RelaxError, RelaxFileError, RelaxNoPdbError, RelaxNoSequenceError, RelaxNoTensorError 38 from lib.io import extract_data, mkdir_nofail, open_write_file, strip, test_binary 39 from pipe_control import angles, diffusion_tensor, pipes, value 40 from pipe_control.interatomic import return_interatom_list 41 from pipe_control.mol_res_spin import exists_mol_res_spin_data, first_residue_num, last_residue_num, residue_loop, return_spin, spin_loop 42 from pipe_control.pipes import check_pipe 43 from pipe_control.spectrometer import loop_frequencies 44 from specific_analyses.model_free.model import determine_model_type 45 46

47 -def __deselect_spins():

48 """Deselect spins with no or too little data, that are overfitting, etc.""" 49 50 # Test if sequence data exists. 51 if not exists_mol_res_spin_data(): 52 raise RelaxNoSequenceError 53 54 # Is structural data required? 55 need_vect = False 56 if hasattr(cdp, 'diff_tensor') and (cdp.diff_tensor.type == 'spheroid' or cdp.diff_tensor.type == 'ellipsoid'): 57 need_vect = True 58 59 # Loop over the sequence. 60 for spin in spin_loop(): 61 # Relaxation data must exist! 62 if not hasattr(spin, 'ri_data'): 63 spin.select = False 64 65 # Require 3 or more relaxation data points. 66 elif len(spin.ri_data) < 3: 67 spin.select = False 68 69 # Require at least as many data points as params to prevent over-fitting. 70 elif hasattr(spin, 'params') and spin.params and len(spin.params) > len(spin.ri_data): 71 spin.select = False

72 73

74 -def create(algor='LM', dir=None, force=False):

75 """Create the Dasha script file 'dasha_script' for controlling the program. 76 77 @keyword algor: The optimisation algorithm to use. This can be the Levenberg-Marquardt algorithm 'LM' or the Newton-Raphson algorithm 'NR'. 78 @type algor: str 79 @keyword dir: The optional directory to place the script into. 80 @type dir: str or None 81 @keyword force: A flag which if True will cause any pre-existing file to be overwritten. 82 @type force: bool 83 """ 84 85 # Test if the current pipe exists. 86 check_pipe() 87 88 # Test if sequence data is loaded. 89 if not exists_mol_res_spin_data(): 90 raise RelaxNoSequenceError 91 92 # Determine the parameter set. 93 model_type = determine_model_type() 94 95 # Test if diffusion tensor data for the data_pipe exists. 96 if model_type != 'local_tm' and not hasattr(cdp, 'diff_tensor'): 97 raise RelaxNoTensorError('diffusion') 98 99 # Test if the PDB file has been loaded (for the spheroid and ellipsoid). 100 if model_type != 'local_tm' and cdp.diff_tensor.type != 'sphere' and not hasattr(cdp, 'structure'): 101 raise RelaxNoPdbError 102 103 # Test the optimisation algorithm. 104 if algor not in ['LM', 'NR']: 105 raise RelaxError("The Dasha optimisation algorithm '%s' is unknown, it should either be 'LM' or 'NR'." % algor) 106 107 # Deselect certain spins. 108 __deselect_spins() 109 110 # Directory creation. 111 if dir == None: 112 dir = pipes.cdp_name() 113 mkdir_nofail(dir, verbosity=0) 114 115 # Calculate the angle alpha of the XH vector in the spheroid diffusion frame. 116 if cdp.diff_tensor.type == 'spheroid': 117 angles.spheroid_frame() 118 119 # Calculate the angles theta and phi of the XH vector in the ellipsoid diffusion frame. 120 elif cdp.diff_tensor.type == 'ellipsoid': 121 angles.ellipsoid_frame() 122 123 # The 'dasha_script' file. 124 script = open_write_file(file_name='dasha_script', dir=dir, force=force) 125 create_script(script, model_type, algor) 126 script.close()

127 128

129 -def create_script(file, model_type, algor):

130 """Create the Dasha script file. 131 132 @param file: The opened file descriptor. 133 @type file: file object 134 @param model_type: The model-free model type. 135 @type model_type: str 136 @param algor: The optimisation algorithm to use. This can be the Levenberg-Marquardt algorithm 'LM' or the Newton-Raphson algorithm 'NR'. 137 @type algor: str 138 """ 139 140 # Delete all data. 141 file.write("# Delete all data.\n") 142 file.write("del 1 10000\n") 143 144 # Nucleus type. 145 file.write("\n# Nucleus type.\n") 146 nucleus = None 147 for spin in spin_loop(): 148 # Skip protons. 149 if spin.isotope == '1H': 150 continue 151 152 # Can only handle one spin type. 153 if nucleus and spin.isotope != nucleus: 154 raise RelaxError("The nuclei '%s' and '%s' do not match, relax can only handle one nucleus type in Dasha." % (nucleus, spin.isotope)) 155 156 # Set the nucleus. 157 if not nucleus: 158 nucleus = spin.isotope 159 160 # Convert the name and write it. 161 if nucleus == '15N': 162 nucleus = 'N15' 163 elif nucleus == '13C': 164 nucleus = 'C13' 165 else: 166 raise RelaxError("Cannot handle the nucleus type '%s' within Dasha." % nucleus) 167 file.write("set nucl %s\n" % nucleus) 168 169 # Number of frequencies. 170 file.write("\n# Number of frequencies.\n") 171 file.write("set n_freq %s\n" % cdp.spectrometer_frq_count) 172 173 # Frequency values. 174 file.write("\n# Frequency values.\n") 175 count = 1 176 for frq in loop_frequencies(): 177 file.write("set H1_freq %s %s\n" % (frq / 1e6, count)) 178 count += 1 179 180 # Set the diffusion tensor. 181 file.write("\n# Set the diffusion tensor.\n") 182 if model_type != 'local_tm': 183 # Sphere. 184 if cdp.diff_tensor.type == 'sphere': 185 file.write("set tr %s\n" % (cdp.diff_tensor.tm / 1e-9)) 186 187 # Spheroid. 188 elif cdp.diff_tensor.type == 'spheroid': 189 file.write('set tr %s\n' % (cdp.diff_tensor.tm / 1e-9)) 190 191 # Ellipsoid. 192 elif cdp.diff_tensor.type == 'ellipsoid': 193 # Get the eigenvales. 194 Dx, Dy, Dz = diffusion_tensor.return_eigenvalues() 195 196 # Geometric parameters. 197 file.write("set tr %s\n" % (cdp.diff_tensor.tm / 1e-9)) 198 file.write("set D1/D3 %s\n" % (Dx / Dz)) 199 file.write("set D2/D3 %s\n" % (Dy / Dz)) 200 201 # Orientational parameters. 202 file.write("set alfa %s\n" % (cdp.diff_tensor.alpha / (2.0 * pi) * 360.0)) 203 file.write("set betta %s\n" % (cdp.diff_tensor.beta / (2.0 * pi) * 360.0)) 204 file.write("set gamma %s\n" % (cdp.diff_tensor.gamma / (2.0 * pi) * 360.0)) 205 206 # Reading the relaxation data. 207 file.write("\n# Reading the relaxation data.\n") 208 file.write("echo Reading the relaxation data.\n") 209 noe_index = 1 210 r1_index = 1 211 r2_index = 1 212 for ri_id in cdp.ri_ids: 213 # NOE. 214 if cdp.ri_type[ri_id] == 'NOE': 215 # Data set number. 216 number = noe_index 217 218 # Data type. 219 data_type = 'noe' 220 221 # Increment the data set index. 222 noe_index = noe_index + 1 223 224 # R1. 225 elif cdp.ri_type[ri_id] == 'R1': 226 # Data set number. 227 number = r1_index 228 229 # Data type. 230 data_type = '1/T1' 231 232 # Increment the data set index. 233 r1_index = r1_index + 1 234 235 # R2. 236 elif cdp.ri_type[ri_id] == 'R2': 237 # Data set number. 238 number = r2_index 239 240 # Data type. 241 data_type = '1/T2' 242 243 # Increment the data set index. 244 r2_index = r2_index + 1 245 246 # Set the data type. 247 if number == 1: 248 file.write("\nread < %s\n" % data_type) 249 else: 250 file.write("\nread < %s %s\n" % (data_type, number)) 251 252 # The relaxation data. 253 for residue in residue_loop(): 254 # Alias the spin. 255 spin = residue.spin[0] 256 257 # Skip deselected spins. 258 if not spin.select: 259 continue 260 261 # Skip and deselect spins for which relaxation data is missing. 262 if len(spin.ri_data) != len(cdp.ri_ids) or spin.ri_data[ri_id] == None: 263 spin.select = False 264 continue 265 266 # Data and errors. 267 file.write("%s %s %s\n" % (residue.num, spin.ri_data[ri_id], spin.ri_data_err[ri_id])) 268 269 # Terminate the reading. 270 file.write("exit\n") 271 272 # Individual residue optimisation. 273 if model_type == 'mf': 274 # Loop over the residues. 275 for residue in residue_loop(): 276 # Alias the spin. 277 spin = residue.spin[0] 278 279 # Skip deselected spins. 280 if not spin.select: 281 continue 282 283 # Get the interatomic data containers. 284 interatoms = return_interatom_list(spin_hash=spin._hash) 285 if len(interatoms) == 0: 286 raise RelaxNoInteratomError 287 elif len(interatoms) > 1: 288 raise RelaxError("Only one interatomic data container, hence dipole-dipole interaction, is supported per spin.") 289 290 # Comment. 291 file.write("\n\n\n# Residue %s\n\n" % residue.num) 292 293 # Echo. 294 file.write("echo Optimisation of residue %s\n" % residue.num) 295 296 # Select the spin. 297 file.write("\n# Select the residue.\n") 298 file.write("set cres %s\n" % residue.num) 299 300 # The angle alpha of the XH vector in the spheroid diffusion frame. 301 if cdp.diff_tensor.type == 'spheroid': 302 file.write("set teta %s\n" % spin.alpha) 303 304 # The angles theta and phi of the XH vector in the ellipsoid diffusion frame. 305 elif cdp.diff_tensor.type == 'ellipsoid': 306 file.write("\n# Setting the spherical angles of the XH vector in the ellipsoid diffusion frame.\n") 307 file.write("set teta %s\n" % spin.theta) 308 file.write("set fi %s\n" % spin.phi) 309 310 # The 'jmode'. 311 if 'ts' in spin.params: 312 jmode = 3 313 elif 'te' in spin.params: 314 jmode = 2 315 elif 's2' in spin.params: 316 jmode = 1 317 318 # Chemical exchange. 319 if 'rex' in spin.params: 320 exch = True 321 else: 322 exch = False 323 324 # Anisotropic diffusion. 325 if cdp.diff_tensor.type == 'sphere': 326 anis = False 327 else: 328 anis = True 329 330 # Axial symmetry. 331 if cdp.diff_tensor.type == 'spheroid': 332 sym = True 333 else: 334 sym = False 335 336 # Set the jmode. 337 file.write("\n# Set the jmode.\n") 338 file.write("set def jmode %s" % jmode) 339 if exch: 340 file.write(" exch") 341 if anis: 342 file.write(" anis") 343 if sym: 344 file.write(" sym") 345 file.write("\n") 346 347 # Parameter default values. 348 file.write("\n# Parameter default values.\n") 349 file.write("reset jmode %s\n" % residue.num) 350 351 # Bond length. 352 file.write("\n# Bond length.\n") 353 file.write("set r_hx %s\n" % (interatoms[0].r / 1e-10)) 354 355 # CSA value. 356 file.write("\n# CSA value.\n") 357 file.write("set csa %s\n" % (spin.csa / 1e-6)) 358 359 # Fix the tf parameter if it isn't in the model. 360 if not 'tf' in spin.params and jmode == 3: 361 file.write("\n# Fix the tf parameter.\n") 362 file.write("fix tf 0\n") 363 364 # Optimisation of all residues. 365 file.write("\n\n\n# Optimisation of all residues.\n") 366 if algor == 'LM': 367 file.write("lmin %s %s" % (first_residue_num(), last_residue_num())) 368 elif algor == 'NR': 369 file.write("min %s %s" % (first_residue_num(), last_residue_num())) 370 371 # Show the results. 372 file.write("\n# Show the results.\n") 373 file.write("echo\n") 374 file.write("show all\n") 375 376 # Write the results. 377 file.write("\n# Write the results.\n") 378 file.write("write s2.out S\n") 379 file.write("write s2f.out Sf\n") 380 file.write("write s2s.out Ss\n") 381 file.write("write te.out te\n") 382 file.write("write tf.out tf\n") 383 file.write("write ts.out ts\n") 384 file.write("write rex.out rex\n") 385 file.write("write chi2.out F\n") 386 387 else: 388 raise RelaxError("Optimisation of the parameter set '%s' currently not supported." % model_type)

389 390

391 -def execute(dir, force, binary):

392 """Execute Dasha. 393 394 @param dir: The optional directory where the script is located. 395 @type dir: str or None 396 @param force: A flag which if True will cause any pre-existing files to be overwritten by Dasha. 397 @type force: bool 398 @param binary: The name of the Dasha binary file. This can include the path to the binary. 399 @type binary: str 400 """ 401 402 # Test the binary file string corresponds to a valid executable. 403 test_binary(binary) 404 405 # The current directory. 406 orig_dir = getcwd() 407 408 # The directory. 409 if dir == None: 410 dir = pipes.cdp_name() 411 if not access(dir, F_OK): 412 raise RelaxDirError('Dasha', dir) 413 414 # Change to this directory. 415 chdir(dir) 416 417 # Catch failures and return to the correct directory. 418 try: 419 # Test if the 'dasha_script' script file exists. 420 if not access('dasha_script', F_OK): 421 raise RelaxFileError('dasha script', 'dasha_script') 422 423 # Python 2.3 and earlier. 424 if Popen == None: 425 raise RelaxError("The subprocess module is not available in this version of Python.") 426 427 # Execute Dasha. 428 pipe = Popen(binary, shell=True, stdin=PIPE, stdout=PIPE, stderr=PIPE, close_fds=False) 429 430 # Get the contents of the script and pump it into Dasha. 431 script = open('dasha_script') 432 lines = script.readlines() 433 script.close() 434 for line in lines: 435 # Encode to a Python 3 byte array. 436 if hasattr(line, 'encode'): 437 line = line.encode() 438 439 # Write out. 440 pipe.stdin.write(line) 441 442 # Close the pipe. 443 pipe.stdin.close() 444 445 # Write to stdout. 446 for line in pipe.stdout.readlines(): 447 # Decode Python 3 byte arrays. 448 if hasattr(line, 'decode'): 449 line = line.decode() 450 451 # Write. 452 sys.stdout.write(line) 453 454 # Write to stderr. 455 for line in pipe.stderr.readlines(): 456 # Decode Python 3 byte arrays. 457 if hasattr(line, 'decode'): 458 line = line.decode() 459 460 # Write. 461 sys.stderr.write(line) 462 463 # Failure. 464 except: 465 # Change back to the original directory. 466 chdir(orig_dir) 467 468 # Reraise the error. 469 raise 470 471 # Change back to the original directory. 472 chdir(orig_dir) 473 474 # Print some blank lines (aesthetics) 475 sys.stdout.write("\n\n")

476 477

478 -def extract(dir):

479 """Extract the data from the Dasha results files. 480 481 @param dir: The optional directory where the results file is located. 482 @type dir: str or None 483 """ 484 485 # Test if sequence data is loaded. 486 if not exists_mol_res_spin_data(): 487 raise RelaxNoSequenceError 488 489 # The directory. 490 if dir == None: 491 dir = pipes.cdp_name() 492 if not access(dir, F_OK): 493 raise RelaxDirError('Dasha', dir) 494 495 # Loop over the parameters. 496 for param in ['s2', 's2f', 's2s', 'te', 'tf', 'ts', 'rex']: 497 # The file name. 498 file_name = dir + sep + param + '.out' 499 500 # Test if the file exists. 501 if not access(file_name, F_OK): 502 raise RelaxFileError('Dasha', file_name) 503 504 # Scaling. 505 if param in ['te', 'tf', 'ts']: 506 scaling = 1e-9 507 elif param == 'rex': 508 scaling = 1.0 / (2.0 * pi * cdp.spectrometer_frq[cdp.ri_ids[0]]) ** 2 509 else: 510 scaling = 1.0 511 512 # Read the values. 513 data = read_results(file=file_name, scaling=scaling) 514 515 # Set the values. 516 for i in range(len(data)): 517 value.set(val=data[i][1], param=param, spin_id=data[i][0]) 518 value.set(val=data[i][0], param=param, spin_id=data[i][0], error=True) 519 520 # Clean up of non-existent parameters (set the parameter to None!). 521 for spin in spin_loop(): 522 # Skip the spin (don't set the parameter to None) if the parameter exists in the model. 523 if param in spin.params: 524 continue 525 526 # Set the parameter to None. 527 setattr(spin, param.lower(), None) 528 529 # Extract the chi-squared values. 530 file_name = dir + sep+'chi2.out' 531 532 # Test if the file exists. 533 if not access(file_name, F_OK): 534 raise RelaxFileError('Dasha', file_name) 535 536 # Read the values. 537 data = read_results(file=file_name) 538 539 # Set the values. 540 for i in range(len(data)): 541 spin = return_spin(spin_id=data[i][0]) 542 spin.chi2 = data[i][1]

543 544

545 -def read_results(file=None, dir=None, scaling=1.0):

546 """Extract the data from the Dasha results file. 547 548 @keyword file: The name of the file to open. 549 @type file: str 550 @keyword dir: The directory containing the file (defaults to the current directory if None). 551 @type dir: str or None 552 @keyword scaling: The parameter scaling factor. 553 @type scaling: float 554 """ 555 556 # Extract the data. 557 data = extract_data(file=file, dir=dir) 558 559 # Remove comments. 560 data = strip(data) 561 562 # Repackage the data as a list of lists of spin ID, value, error. 563 new_data = [] 564 for i in range(len(data)): 565 spin_id = ':%s@N' % data[i][0] 566 value = float(data[i][1]) * scaling 567 error = float(data[i][2]) * scaling 568 new_data.append([spin_id, value, error]) 569 570 # Return the data. 571 return new_data

572

Source Code for Module pipe_control.dasha