Coverage for biobb_analysis/ambertools/common.py: 63%
274 statements
« prev ^ index » next coverage.py v7.15.3, created at 2026-08-07 12:30 +0000
1""" Common functions for package biobb_analysis.ambertools """
2from pathlib import Path, PurePath
3import zipfile
4import shutil
5from biobb_common.tools import file_utils as fu
8def check_top_path(path, out_log, classname):
9 """ Checks topology input file """
10 orig_path = path
11 if not Path(path).exists():
12 fu.log(classname + ': Unexisting topology input file, exiting', out_log)
13 raise SystemExit(classname + ': Unexisting topology input file')
14 file_extension = PurePath(path).suffix
15 if not is_valid_topology(file_extension[1:]):
16 fu.log(classname + ': Format %s in topology input file is not compatible' % file_extension[1:], out_log)
17 raise SystemExit(classname + ': Format %s in topology input file is not compatible' % file_extension[1:])
18 if zipfile.is_zipfile(path):
19 top_file = fu.unzip_top(zip_file=path, out_log=out_log)
20 path = top_file
21 return path, orig_path
24def check_traj_path(path, out_log, classname):
25 """ Checks trajectory input file """
26 if not Path(path).exists():
27 fu.log(classname + ': Unexisting trajectory input file, exiting', out_log)
28 raise SystemExit(classname + ': Unexisting trajectory input file')
29 file_extension = PurePath(path).suffix
30 if not is_valid_trajectory(file_extension[1:]):
31 fu.log(classname + ': Format %s in trajectory input file is not compatible' % file_extension[1:], out_log)
32 raise SystemExit(classname + ': Format %s in trajectory input file is not compatible' % file_extension[1:])
33 return path
36def check_out_path(path, out_log, classname):
37 """ Checks if output folder exists """
38 if PurePath(path).parent and not Path(PurePath(path).parent).exists():
39 fu.log(classname + ': Unexisting output folder, exiting', out_log)
40 raise SystemExit(classname + ': Unexisting output folder')
41 return path
44def get_parameters(properties, type, classname, out_log):
45 """ Gets in_parameters and out_parameters """
46 if not properties.get(type, dict()):
47 fu.log('No %s parameters provided' % type, out_log)
48 return get_default_value(classname)[type]
50 return {k: v for k, v in properties.get(type, dict()).items()}
53def get_binary_path(properties, type):
54 """ Gets binary path """
55 return properties.get(type, get_default_value(type))
58def check_in_path(path, out_log, classname):
59 """ Checks input instructions file """
60 if not Path(path).exists():
61 fu.log(classname + ': Unexisting input instructions file, exiting', out_log)
62 raise SystemExit(classname + ': Unexisting input instructions file')
64 # check syntax for instructions file
65 syntax_instructions = True
66 err_instructions = ''
67 if 'parm ' not in open(path).read():
68 syntax_instructions = False
69 err_instructions += 'No topology provided'
70 if 'trajin ' not in open(path).read():
71 syntax_instructions = False
72 err_instructions += ' No input trajectory provided'
73 if not syntax_instructions:
74 fu.log(classname + ': Incorrect syntax for instructions file: %s, exiting' % err_instructions, out_log)
75 raise SystemExit(classname + ': Incorrect syntax for instructions file: %s, exiting' % err_instructions)
78def get_default_value(key):
79 """ Gives default values according to the given key """
80 default_values = {
81 "start": 1,
82 "end": -1,
83 "step": 1,
84 "snapshot": 1,
85 "format": "netcdf",
86 "mask": "all-atoms",
87 "reference": "first",
88 "average": "MyAvg",
89 "instructions_file": "instructions.in",
90 "binary_path": "cpptraj",
91 "algorithm": "hieragglo",
92 "metric": "rms",
93 "linkage": "linkage",
94 "clusters": 10,
95 "minpoints": 4,
96 "cluster_set": "MyClusters",
97 # default conf for Average
98 "Average": {
99 "in_parameters": {
100 "start": 1,
101 "end": -1,
102 "step": 1,
103 "mask": "all-atoms"
104 },
105 "out_parameters": {
106 "format": "pdb"
107 }
108 },
109 # default conf for Bfactor
110 "Bfactor": {
111 "in_parameters": {
112 "start": 1,
113 "end": -1,
114 "step": 1,
115 "mask": "all-atoms",
116 "reference": "first"
117 }
118 },
119 # default conf for Cluster
120 "Cluster": {
121 "in_parameters": {
122 "start": 1,
123 "end": -1,
124 "step": 1,
125 "mask": "all-atoms"
126 }
127 },
128 # default conf for Convert
129 "Convert": {
130 "in_parameters": {
131 "start": 1,
132 "end": -1,
133 "step": 1,
134 "mask": "all-atoms"
135 },
136 "out_parameters": {
137 "format": "netcdf"
138 }
139 },
140 # default conf for Dry
141 "Dry": {
142 "in_parameters": {
143 "start": 1,
144 "end": -1,
145 "step": 1,
146 "mask": "all-atoms"
147 },
148 "out_parameters": {
149 "format": "netcdf"
150 }
151 },
152 # default conf for Image
153 "Image": {
154 "in_parameters": {
155 "start": 1,
156 "end": -1,
157 "step": 1,
158 "mask": "all-atoms"
159 },
160 "out_parameters": {
161 "format": "netcdf"
162 }
163 },
164 # default conf for Mask
165 "Mask": {
166 "in_parameters": {
167 "start": 1,
168 "end": -1,
169 "step": 1,
170 "mask": "all-atoms"
171 },
172 "out_parameters": {
173 "format": "netcdf"
174 }
175 },
176 # default conf for Rgyr
177 "Rgyr": {
178 "in_parameters": {
179 "start": 1,
180 "end": -1,
181 "step": 1,
182 "mask": "all-atoms"
183 }
184 },
185 # default conf for Rms
186 "Rms": {
187 "in_parameters": {
188 "start": 1,
189 "end": -1,
190 "step": 1,
191 "mask": "all-atoms",
192 "reference": "first"
193 }
194 },
195 # default conf for Rmsf
196 "Rmsf": {
197 "in_parameters": {
198 "start": 1,
199 "end": -1,
200 "step": 1,
201 "mask": "all-atoms",
202 "reference": "first"
203 }
204 },
205 # default conf for Slice
206 "Slice": {
207 "in_parameters": {
208 "start": 1,
209 "end": -1,
210 "step": 1,
211 "mask": "all-atoms"
212 },
213 "out_parameters": {
214 "format": "netcdf"
215 }
216 },
217 # default conf for Snapshot
218 "Snapshot": {
219 "in_parameters": {
220 "snapshot": 12,
221 "mask": "all-atoms"
222 },
223 "out_parameters": {
224 "format": "pdb"
225 }
226 },
227 # default conf for Strip
228 "Strip": {
229 "in_parameters": {
230 "start": 1,
231 "end": -1,
232 "step": 1,
233 "mask": "all-atoms"
234 },
235 "out_parameters": {
236 "format": "netcdf"
237 }
238 }
239 }
241 return default_values[key]
244def is_valid_topology(ext):
245 """ Checks if trajectory format is compatible with Cpptraj """
246 formats = 'top', 'pdb', 'prmtop', 'parmtop', 'zip'
247 return ext in formats
250def is_valid_trajectory(traj):
251 """ Checks if trajectory format is compatible with Cpptraj """
252 formats = 'mdcrd', 'crd', 'cdf', 'netcdf', 'nc', 'restart', 'ncrestart', 'restartnc', 'dcd', 'charmm', 'cor', 'pdb', 'mol2', 'trr', 'gro', 'binpos', 'xtc', 'cif', 'arc', 'sqm', 'sdf', 'conflib'
253 return traj in formats
256def is_valid_reference(ref):
257 """ Checks if reference is correct """
258 references = 'first', 'average', 'experimental'
259 return ref in references
262def is_valid_cluster_algorithm(algorithm):
263 """ Checks if clustering algorithm is compatible with Cpptraj """
264 algorithms = 'hieragglo', 'dbscan', 'kmeans', 'dpeaks'
265 return algorithm in algorithms
268def is_valid_cluster_metric(metric):
269 """ Checks if clustering distance metric is compatible with Cpptraj """
270 metrics = 'rms', 'srmsd', 'dme', 'qrmsd'
271 return metric in metrics
274def is_valid_cluster_linkage(linkage):
275 """ Checks if hierarchical agglomerative linkage is compatible with Cpptraj """
276 linkages = 'linkage', 'averagelinkage', 'complete'
277 return linkage in linkages
280def get_cluster_algorithm(algorithm, clusters, epsilon, minpoints, linkage, out_log):
281 """ Return string with the clustering algorithm keywords """
282 if not algorithm or algorithm == 'None':
283 algorithm = get_default_value('algorithm')
284 fu.log('No algorithm provided in configuration file, assigned default value: %s' % get_default_value('algorithm'), out_log)
286 if not is_valid_cluster_algorithm(algorithm):
287 fu.log('Algorithm %s is not compatible, assigned default value: %s' % (algorithm, get_default_value('algorithm')), out_log)
288 algorithm = get_default_value('algorithm')
290 if algorithm == 'hieragglo':
291 if not is_valid_cluster_linkage(linkage):
292 fu.log('Linkage %s is not compatible, assigned default value: %s' % (linkage, get_default_value('linkage')), out_log)
293 linkage = get_default_value('linkage')
294 keywords = ['hieragglo', 'clusters', str(clusters), linkage]
295 # epsilon is an additional stop condition for hieragglo, only added when provided
296 if epsilon is not None and epsilon > 0:
297 keywords += ['epsilon', str(epsilon)]
298 elif algorithm == 'kmeans':
299 keywords = ['kmeans', 'clusters', str(clusters)]
300 else:
301 # both dbscan and dpeaks need epsilon
302 if epsilon is None or epsilon <= 0:
303 fu.log('Algorithm %s needs a positive epsilon, assigned default value: 1.0' % algorithm, out_log)
304 epsilon = 1.0
305 if algorithm == 'dbscan':
306 keywords = ['dbscan', 'minpoints', str(minpoints), 'epsilon', str(epsilon)]
307 else:
308 keywords = ['dpeaks', 'epsilon', str(epsilon), 'choosepoints', 'auto']
310 return ' '.join(keywords)
313def get_cluster_metric(metric, mass, nofit, metric_mask, out_log):
314 """ Return string with the clustering distance metric keywords """
315 if not metric or metric == 'None':
316 metric = get_default_value('metric')
317 fu.log('No metric provided in configuration file, assigned default value: %s' % get_default_value('metric'), out_log)
319 if not is_valid_cluster_metric(metric):
320 fu.log('Metric %s is not compatible, assigned default value: %s' % (metric, get_default_value('metric')), out_log)
321 metric = get_default_value('metric')
323 keywords = [metric]
324 if mass:
325 keywords.append('mass')
326 if nofit:
327 keywords.append('nofit')
329 if not metric_mask or metric_mask == 'None':
330 metric_mask = get_default_value('mask')
331 fu.log('No metric_mask provided in configuration file, assigned default value: %s' % get_default_value('mask'), out_log)
332 # the metric mask selects atoms, it does not remove them from the clustered coordinates
333 keywords.append(get_mask(metric_mask, out_log))
335 return ' '.join(keywords)
338def get_traj_format(path, out_log):
339 """ Return the Cpptraj trajectory format matching the extension of the given path """
340 format = PurePath(path).suffix[1:]
341 if not is_valid_trajectory(format):
342 fu.log('Format %s is not compatible, assigned default value: %s' % (format, get_default_value('format')), out_log)
343 format = get_default_value('format')
345 return format
348def get_mask_atoms(key):
349 """ Gives mask atoms according to the given key """
350 masks = {
351 "c-alpha": "@CA",
352 "backbone": "@C,CA,N,O,C3',O3',C4',C5',O5',P",
353 "all-atoms": ":*",
354 "heavy-atoms": "!@H*,1H*,2H*,3H*",
355 "side-chain": "!@CA,C,N,O,H,HA,C3',O3',C4',C5',O5',P",
356 "solute": "!:WAT,HOH,SOL,TIP3,TP3,SOD,CLA,Na+,Cl-,NA,CL,K+,K",
357 "ions": ":SOD,CLA,Na+,Cl-,NA,CL,K+,K",
358 "solvent": ":WAT,HOH,SOL,TIP3,TP3"
359 }
361 # if key incorrect, return default value and message
362 if key in masks:
363 return masks[key], None
364 else:
365 return key, None # Allow for Amber mask
368def get_in_parameters(list, out_log, type='None'):
369 """ Return string with input parameters """
370 # if strip or mask, no mandatory trajin parameters
371 if type == 'strip' or type == 'mask':
372 start = '' if 'start' not in list else str(list['start'])
373 end = '' if 'end' not in list else str(list['end'])
374 step = '' if 'step' not in list else str(list['step'])
375 if (not start or start == 'None') and (not end or end == 'None') and (not step or step == 'None'):
376 return ''
377 else:
378 if not start:
379 start = str(get_default_value("start"))
380 fu.log('No start value provided in configuration file or incorrect format, assigned default value: %s' % get_default_value('start'), out_log)
381 if not end:
382 end = str(get_default_value("end"))
383 fu.log('No end value provided in configuration file or incorrect format, assigned default value: %s' % get_default_value('end'), out_log)
384 if not step:
385 step = str(get_default_value("step"))
386 fu.log('No step value provided in configuration file or incorrect format, assigned default value: %s' % get_default_value('step'), out_log)
387 else:
388 # check if trajin parameters are provided and have correct format
389 if (type == 'snapshot'):
390 snapshot = str(get_default_value("snapshot")) if 'snapshot' not in list else str(list['snapshot'])
391 if 'snapshot' not in list or not list['snapshot'] or not isinstance(list['snapshot'], int):
392 snapshot = str(get_default_value("snapshot"))
393 fu.log('No snapshot value provided in configuration file or incorrect format, assigned default value: %s' % get_default_value('snapshot'), out_log)
394 start = snapshot
395 end = snapshot
396 step = '1'
397 else:
398 start = str(get_default_value("start")) if 'start' not in list else str(list['start'])
399 if 'start' not in list or not list['start'] or not isinstance(list['start'], int):
400 start = str(get_default_value("start"))
401 fu.log('No start value provided in configuration file or incorrect format, assigned default value: %s' % get_default_value('start'), out_log)
402 end = str(get_default_value("end")) if 'end' not in list else str(list['end'])
403 if 'end' not in list or not list['end'] or not isinstance(list['end'], int):
404 end = str(get_default_value("end"))
405 fu.log('No end value provided in configuration file or incorrect format, assigned default value: %s' % get_default_value('end'), out_log)
406 step = str(get_default_value("step")) if 'step' not in list else str(list['step'])
407 if 'step' not in list or not list['step'] or not isinstance(list['step'], int):
408 step = str(get_default_value("step"))
409 fu.log('No step value provided in configuration file or incorrect format, assigned default value: %s' % get_default_value('step'), out_log)
411 # checking start <= end
412 if end != '-1' and start > end:
413 fu.log('End must be -1 (indicating the end of the trajectory) or more or equal than start. Your values are start: %s, end: %s' % (start, end), out_log)
414 raise SystemExit('End must be -1 (indicating the end of the trajectory) or greater or equal than start. Your values are start: %s, end: %s' % (start, end))
416 return start + " " + end + " " + step
419def setup_structure(out_log):
420 """ Sets up the structure """
421 instructions_list = []
422 mask_atoms = get_mask('heavy-atoms', out_log)
423 instructions_list.append('center ' + mask_atoms + ' origin')
424 instructions_list.append('autoimage')
425 instructions_list.append('rms first ' + mask_atoms)
426 mask_solvent = get_mask('solvent', out_log)
427 mask_ions = get_mask('ions', out_log)
428 instructions_list.append('strip ' + mask_solvent + ',' + mask_ions[1:])
430 return instructions_list
433def get_negative_mask(key, out_log):
434 """ Gives the negative mask according to the given key """
435 atoms, msg = get_mask_atoms(key)
436 if atoms[0] == '!':
437 mask = atoms[1:]
438 else:
439 mask = '!' + atoms
440 # if mask incorrect, give message
441 if msg:
442 fu.log(msg, out_log)
444 return mask
447def get_mask(key, out_log):
448 """ Gives mask according to the given key """
449 mask, msg = get_mask_atoms(key)
450 # if mask incorrect, give message
451 if msg:
452 fu.log(msg, out_log)
454 return mask
457def get_reference(ref, output_cpptraj_path, input_exp_path, mask, output, classname, out_log):
458 """ Gives reference instructions according to the given key """
459 instructions_list = []
460 if not ref or ref == 'None':
461 ref = get_default_value('reference')
462 fu.log('No reference provided in configuration file, assigned default value: %s' % get_default_value('reference'), out_log)
464 if not is_valid_reference(ref):
465 fu.log('Reference %s is not compatible, assigned default value: %s' % (ref, get_default_value('reference')), out_log)
466 ref = get_default_value('reference')
468 if ref == 'first':
469 if output:
470 instructions_list.append('rms first out ' + output_cpptraj_path)
471 else:
472 instructions_list.append('rms first')
474 if ref == 'average':
475 instructions_list.append('average crdset ' + get_default_value('average'))
476 instructions_list.append('run')
477 if output:
478 instructions_list.append('rms ref ' + get_default_value('average') + ' ' + mask + ' out ' + output_cpptraj_path)
479 else:
480 instructions_list.append('rms ref ' + get_default_value('average') + ' ' + mask)
482 if ref == 'experimental':
483 if not input_exp_path:
484 fu.log('No experimental structure provided, exiting', out_log)
485 raise SystemExit(classname + ': input_exp_path is mandatory')
486 instructions_list.append('parm ' + input_exp_path + ' noconect [exp]')
487 solute, msg = get_mask_atoms('solute')
488 instructions_list.append('reference ' + input_exp_path + ' ' + solute + ' parm [exp]')
489 backbone, msg = get_mask_atoms('backbone')
490 if output:
491 instructions_list.append('rms reference ' + mask + ' out ' + output_cpptraj_path)
492 else:
493 instructions_list.append('rms reference ' + mask)
495 return instructions_list
498def get_reference_rms(ref, output_cpptraj_path, input_exp_path, mask, output, classname, out_log, nofit=False, norotate=False, nomod=False):
499 """ Gives reference instructions according to the given key """
500 instructions_list = []
501 if not ref or ref == 'None':
502 ref = get_default_value('reference')
503 fu.log('No reference provided in configuration file, assigned default value: %s' % get_default_value('reference'), out_log)
505 if not is_valid_reference(ref):
506 fu.log('Reference %s is not compatible, assigned default value: %s' % (ref, get_default_value('reference')), out_log)
507 ref = get_default_value('reference')
509 flags = []
510 if nofit:
511 flags.append("nofit")
512 if norotate:
513 flags.append("norotate")
514 if nomod:
515 flags.append("nomod")
517 flags_str = " ".join(flags)
519 if ref == 'first':
520 if output:
521 instructions_list.append('rms first out ' + output_cpptraj_path + f' {flags_str}')
522 else:
523 instructions_list.append('rms first' + f' {flags_str}')
525 if ref == 'average':
526 instructions_list.append('average crdset ' + get_default_value('average'))
527 instructions_list.append('run')
528 if output:
529 instructions_list.append('rms ref ' + get_default_value('average') + ' ' + mask + ' out ' + output_cpptraj_path + f' {flags_str}')
530 else:
531 instructions_list.append('rms ref ' + get_default_value('average') + ' ' + mask + f' {flags_str}')
533 if ref == 'experimental':
534 if not input_exp_path:
535 fu.log('No experimental structure provided, exiting', out_log)
536 raise SystemExit(classname + ': input_exp_path is mandatory')
537 instructions_list.append('parm ' + input_exp_path + ' noconect [exp]')
538 solute, msg = get_mask_atoms('solute')
539 instructions_list.append('reference ' + input_exp_path + ' ' + solute + ' parm [exp]')
540 backbone, msg = get_mask_atoms('backbone')
541 if output:
542 instructions_list.append('rms reference ' + mask + ' out ' + output_cpptraj_path + f' {flags_str}')
543 else:
544 instructions_list.append('rms reference ' + mask + f' {flags_str}')
546 return instructions_list
549def get_out_parameters(list, out_log):
550 """ Return string with output parameters """
551 format = list['format']
552 # msg = None
553 # check if format provided
554 if not format:
555 format = get_default_value('format')
556 fu.log('No format provided in configuration file, assigned default value: %s' % get_default_value('format'), out_log)
557 # check if valid format
558 if not is_valid_trajectory(format):
559 fu.log('Format %s is not compatible, assigned default value: %s' % (format, get_default_value('format')), out_log)
560 format = get_default_value('format')
562 return format
565def copy_instructions_file_to_container(instructions_file, unique_dir):
566 shutil.copy2(instructions_file, unique_dir)
569def remove_tmp_files(list, remove_tmp, out_log, input_top_path_orig=None, input_top_path=None):
570 """ Removes temporal files generated by the wrapper """
571 tmp_files = list
572 if zipfile.is_zipfile(str(input_top_path_orig)):
573 tmp_files.append(PurePath(str(input_top_path)).parent)
575 if remove_tmp:
576 removed_files = [f for f in tmp_files if fu.rm(f)]
577 fu.log('Removed: %s' % str(removed_files), out_log)