lib.spectrum.nmrpipe

1 ############################################################################### 2 # # 3 # Copyright (C) 2013-2014 Troels E. Linnet # 4 # Copyright (C) 2014 Edward d'Auvergne # 5 # # 6 # This file is part of the program relax (http://www.nmr-relax.com). # 7 # # 8 # This program is free software: you can redistribute it and/or modify # 9 # it under the terms of the GNU General Public License as published by # 10 # the Free Software Foundation, either version 3 of the License, or # 11 # (at your option) any later version. # 12 # # 13 # This program is distributed in the hope that it will be useful, # 14 # but WITHOUT ANY WARRANTY; without even the implied warranty of # 15 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # 16 # GNU General Public License for more details. # 17 # # 18 # You should have received a copy of the GNU General Public License # 19 # along with this program. If not, see <http://www.gnu.org/licenses/>. # 20 # # 21 ############################################################################### 22 23 # Module docstring. 24 """Module containing functions for handling NMRPipe SeriesTab files.""" 25 26 27 # Python module imports. 28 import re 29 from glob import glob 30 from os import sep 31 from os.path import abspath 32 subprocess_module = False 33 try: 34 import subprocess 35 subprocess_module = True 36 except ImportError: 37 pass 38 from warnings import warn 39 40 # relax module imports. 41 from lib.errors import RelaxError 42 from lib.io import file_root, get_file_path, open_write_file, sort_filenames, write_data 43 from lib.warnings import RelaxWarning 44 45

46 -def read_seriestab(peak_list=None, file_data=None, int_col=None):

47 """Extract the intensity information from the NMRPipe SeriesTab peak intensity file. 48 49 @keyword peak_list: The peak list object to place all data into. 50 @type peak_list: lib.spectrum.objects.Peak_list instance 51 @keyword file_data: The data extracted from the file converted into a list of lists. 52 @type file_data: list of lists of str 53 @keyword int_col: The column which to multiply the peak intensity data (used by the SeriesTab intensity file format). 54 @type int_col: int 55 @raises RelaxError: When the expected peak intensity is not a float. 56 """ 57 58 # Set start variables. 59 modeline = False 60 mode = False 61 varsline = False 62 header = False 63 64 # Loop over lines, to extract variables and find header size. 65 line_nr = 0 66 for line in file_data: 67 if len(line) > 0: 68 if line[0] == 'REMARK' and line[1] == 'Mode:': 69 modeline = line[2:] 70 mode = modeline[0] 71 elif line[0] == 'VARS': 72 varsline = line[1:] 73 elif line[0] == '1': 74 header = line_nr 75 break 76 line_nr += 1 77 78 # Raise RelaxError, if the MODE is not found. 79 if not (modeline and mode): 80 raise RelaxError("MODE not detected. Expecting line 2:\nREMARK Mode: Summation") 81 82 # Raise RelaxError, if the VARS line is not found. 83 if not (varsline): 84 raise RelaxError("VARS not detected. Expecting line 8:\nVARS INDEX X_AXIS Y_AXIS X_PPM Y_PPM VOL ASS Z_A0") 85 86 # Raise RelaxError, if the header size is not found. 87 if not header: 88 raise RelaxError("'1' not detected in start of line. Cannot determine header size.") 89 90 # Find index of assignment ASS. 91 ass_i = varsline.index('ASS') 92 93 # Chemical shifts preparation. 94 w1_col = None 95 w2_col = None 96 97 # Find index of chemical shift Y_PPM which in sparky is w1. 98 w1_col = varsline.index('Y_PPM') 99 100 # Find index of chemical shift X_PPM which in sparky is w2. 101 w2_col = varsline.index('X_PPM') 102 103 # Make a regular search for Z_A entries. 104 Z_A = re.compile("Z_A*") 105 spectra = list(filter(Z_A.search, varsline)) 106 107 # Find index of Z_A entries. 108 spectra_i = [] 109 for y in spectra: 110 spectra_i.append(varsline.index(y)) 111 112 # Remove the header. 113 file_data = file_data[header:] 114 115 # Loop over the file data. 116 for line in file_data: 117 # Skip non-assigned peaks. 118 if line[ass_i] == '?-?': 119 continue 120 121 # First split by the 2D separator. 122 assign1, assign2 = re.split('-', line[ass_i]) 123 124 # The assignment of the first dimension. 125 row1 = re.split('([a-zA-Z]+)', assign1) 126 name1 = row1[-2] + row1[-1] 127 128 # The assignment of the second dimension. 129 row2 = re.split('([a-zA-Z]+)', assign2) 130 name2 = row2[-2] + row2[-1] 131 132 # Get the residue number for dimension 1. 133 got_res_num1 = True 134 try: 135 res_num1 = int(row1[-3]) 136 except: 137 got_res_num1 = False 138 raise RelaxError("Improperly formatted NMRPipe SeriesTab file, cannot process the residue number for dimension 1 in assignment: %s." % line[0]) 139 140 # Get the residue number for dimension 2. 141 try: 142 res_num2 = int(row2[-3]) 143 except: 144 # We cannot always expect dimension 2 to have residue number. 145 if got_res_num1: 146 res_num2 = res_num1 147 else: 148 res_num2 = None 149 warn(RelaxWarning("Improperly formatted NMRPipe SeriesTab file, cannot process the residue number for dimension 2 in assignment: %s. Setting residue number to %s." % (line[0], res_num2))) 150 151 # The residue name for dimension 1. 152 got_res_name1 = True 153 try: 154 res_name1 = row1[-4] 155 except: 156 got_res_name1 = False 157 res_name1 = None 158 warn(RelaxWarning("Improperly formatted NMRPipe SeriesTab file, cannot process the residue name for dimension 1 in assignment: %s. Setting residue name to %s." % (line[0], res_name1))) 159 160 # The residue name for dimension 2. 161 try: 162 res_name2 = row2[-4] 163 except: 164 # We cannot always expect dimension 2 to have residue name. 165 if got_res_name1: 166 res_name2 = res_name1 167 else: 168 res_name2 = None 169 warn(RelaxWarning("Improperly formatted NMRPipe SeriesTab file, cannot process the residue name for dimension 2 in assignment: %s. Setting residue name to %s." % (line[0], res_name2))) 170 171 # Get the intensities. 172 try: 173 # Loop over the spectra. 174 intensities = [] 175 for i in range(len(spectra)): 176 # The intensity is given by column multiplication. 177 intensities.append(float(line[spectra_i[i]])*float(line[5])) 178 179 # Bad data. 180 except ValueError: 181 raise RelaxError("The peak intensity value %s from the line %s is invalid." % (intensity, line)) 182 183 # Chemical shifts. 184 w1 = None 185 w2 = None 186 if w1_col != None: 187 try: 188 w1 = float(line[w1_col]) 189 except ValueError: 190 raise RelaxError("The chemical shift from the line %s is invalid." % line) 191 if w2_col != None: 192 try: 193 w2 = float(line[w2_col]) 194 except ValueError: 195 raise RelaxError("The chemical shift from the line %s is invalid." % line) 196 197 # Add the assignment to the peak list object. 198 peak_list.add(res_nums=[res_num1, res_num2], res_names=[res_name1, res_name2], spin_names=[name1, name2], shifts=[w1, w2], intensity=intensities, intensity_name=spectra)

199 200

201 -def show_apod_extract(file_name=None, dir=None, path_to_command='showApod'):

202 """Extract showApod information for spectrum fourier transformed with NMRPipe. 203 204 @keyword file: The filename of the NMRPipe fourier transformed file. 205 @type file: str 206 @keyword dir: The directory where the file is located. 207 @type dir: str 208 @keyword path_to_command: If showApod not in PATH, then specify absolute path as: /path/to/showApod 209 @type path_to_command: str 210 @return: The output from showApod as list of lines. 211 @rtype: list of lines 212 """ 213 214 # Get the file path. 215 file_path = get_file_path(file_name=file_name, dir=dir) 216 217 if not subprocess_module: 218 raise RelaxError("Python module 'subprocess' not found, cannot call showApod.") 219 220 # Call function. 221 Temp = subprocess.Popen([path_to_command, file_path], stdout=subprocess.PIPE) 222 223 # Communicate with program, and get outout and exitcode. 224 (output, errput) = Temp.communicate() 225 226 # Wait for finish and get return code. 227 return_value = Temp.wait() 228 229 # Python 3 support - convert byte arrays to text. 230 if hasattr(output, 'decode'): 231 output = output.decode() 232 233 return output.splitlines()

234 235

236 -def show_apod_rmsd(file_name=None, dir=None, path_to_command='showApod'):

237 """Extract showApod 'Noise Std Dev' for spectrum fourier transformed with NMRPipe. 238 239 @keyword file: The filename of the NMRPipe fourier transformed file. 240 @type file: str 241 @keyword dir: The directory where the file is located. 242 @type dir: str 243 @keyword path_to_command: If showApod not in PATH, then specify absolute path as: /path/to/showApod 244 @type path_to_command: str 245 @return: The Noise Std Dev from line: 'REMARK Automated Noise Std Dev in Processed Data' 246 @rtype: float 247 """ 248 249 # Call extract function. 250 show_apod_lines = show_apod_extract(file_name=file_name, dir=dir, path_to_command=path_to_command) 251 252 # Loop over the lines 253 found = False 254 for line in show_apod_lines: 255 # Look for line with this remark. 256 if line[:49] == 'REMARK Automated Noise Std Dev in Processed Data:': 257 # The rest of the line is the rmsd. 258 rmsd = float(line[49:].split()[0]) 259 return rmsd 260 261 if not found: 262 print(show_apod_lines) 263 raise RelaxError("Could not find the line: 'REMARK Automated Noise Std Dev in Processed Data:', from the output of showApod.")

264 265

266 -def show_apod_rmsd_to_file(file_name=None, dir=None, path_to_command='showApod', outdir=None, force=False):

267 """Extract showApod 'Noise Std Dev' from showApod, and write to file with same filename and ending '.rmsd' 268 269 @keyword file: The filename of the NMRPipe fourier transformed file. 270 @type file: str 271 @keyword dir: The directory where the file is located. 272 @type dir: str 273 @keyword path_to_command: If showApod not in PATH, then specify absolute path as: /path/to/showApod 274 @type path_to_command: str 275 @keyword outdir: The directory where to write the file. If 'None', then write in same directory. 276 @type outdir: str 277 @param force: Boolean argument which if True causes the file to be overwritten if it already exists. 278 @type force: bool 279 @return: Write the 'Noise Std Dev' from showApod to a file with same file filename, with ending '.rmsd'. This will be a file path. 280 @rtype: str 281 """ 282 283 # Call extract function. 284 apod_rmsd = show_apod_rmsd(file_name=file_name, dir=dir, path_to_command=path_to_command) 285 286 # Get the filename striped of extension details. 287 file_name_root = file_root(file_name) 288 289 # Define extension. 290 extension = ".rmsd" 291 292 # Define file name for writing. 293 file_name_out = file_name_root + extension 294 295 # Define folder to write to. 296 if outdir == None: 297 write_outdir = dir 298 else: 299 write_outdir = outdir 300 301 # Open file for writing, 302 wfile, wfile_path = open_write_file(file_name=file_name_out, dir=write_outdir, force=force, verbosity=1, return_path=True) 303 304 # Write to file. 305 out_write_data = [['%s'%apod_rmsd]] 306 307 # Write data 308 write_data(out=wfile, headings=None, data=out_write_data, sep=None) 309 310 # Close file. 311 wfile.close() 312 313 # Return path to file. 314 return wfile_path

315 316

317 -def show_apod_rmsd_dir_to_files(file_ext='.ft2', dir=None, path_to_command='showApod', outdir=None, force=False):

318 """Searches 'dir' for files with extension 'file_ext'. Extract showApod 'Noise Std Dev' from showApod, and write to file with same filename and ending '.rmsd'. 319 320 @keyword file_ext: The extension for files which is NMRPipe fourier transformed file. 321 @type file_ext: str 322 @keyword dir: The directory where the files is located. 323 @type dir: str 324 @keyword path_to_command: If showApod not in PATH, then specify absolute path as: /path/to/showApod 325 @type path_to_command: str 326 @keyword outdir: The directory where to write the files. If 'None', then write in same directory. 327 @type outdir: str 328 @param force: Boolean argument which if True causes the file to be overwritten if it already exists. 329 @type force: bool 330 @return: Write the 'Noise Std Dev' from showApod to a file with same file filename, with ending '.rmsd'. This will be a list of file paths. 331 @rtype: list of str 332 """ 333 334 # First get correct dir, no matter if dir is specified with or without system folder separator. 335 dir_files = abspath(dir) 336 337 # Define glop pattern. 338 glob_pat = '*%s' % file_ext 339 340 # Get a list of files which math the file extension. 341 file_list = glob(dir_files + sep + glob_pat) 342 343 # Now sort into Alphanumeric order. 344 file_list_sorted = sort_filenames(filenames=file_list, rev=False) 345 346 # Loop over the files. 347 rmsd_files = [] 348 for ft_file in file_list_sorted: 349 # Write rmsd to file, and get file path to file. 350 rmsd_file = show_apod_rmsd_to_file(file_name=ft_file, path_to_command=path_to_command, outdir=outdir, force=force) 351 352 # Collect file path. 353 rmsd_files.append(rmsd_file) 354 355 return rmsd_files

356

Source Code for Module lib.spectrum.nmrpipe