biobb_gromacs/gromacs_extra/ndx2resttop.py source
1 #!/usr/bin/env python3
2
3 """Module containing the Ndx2resttop class and the command line interface."""
4 import fnmatch
5 import argparse
6 from typing import Optional
7 from typing import Any
8 from pathlib import Path
9 from biobb_common.generic.biobb_object import BiobbObject
10 from biobb_common.configuration import settings
11 from biobb_common.tools import file_utils as fu
12 from biobb_common.tools.file_utils import launchlogger
13
14
15 class Ndx2resttop(BiobbObject):
16 """
17 | biobb_gromacs Ndx2resttop
18 | Generate a restrained topology from an index NDX file.
19 | This module automatizes the process of restrained topology generation starting from an index NDX file.
20
21 Args:
22 input_ndx_path (str): Path to the input NDX index file. File type: input. `Sample file <https://github.com/bioexcel/biobb_gromacs/raw/master/biobb_gromacs/test/data/gromacs_extra/ndx2resttop.ndx>`_. Accepted formats: ndx (edam:format_2033).
23 input_top_zip_path (str): Path the input TOP topology in zip format. File type: input. `Sample file <https://github.com/bioexcel/biobb_gromacs/raw/master/biobb_gromacs/test/data/gromacs_extra/ndx2resttop.zip>`_. Accepted formats: zip (edam:format_3987).
24 output_top_zip_path (str): Path the output TOP topology in zip format. File type: output. `Sample file <https://github.com/bioexcel/biobb_gromacs/raw/master/biobb_gromacs/test/reference/gromacs_extra/ref_ndx2resttop.zip>`_. Accepted formats: zip (edam:format_3987).
25 properties (dict - Python dictionary object containing the tool parameters, not input/output files):
26 * **posres_names** (*list*) - ("CUSTOM_POSRES") String with names of the position restraints to be included in the topology file separated by spaces. If provided it should match the length of the ref_rest_chain_triplet_list.
27 * **force_constants** (*str*) - ("500 500 500") Array of three floats defining the force constants.
28 * **ref_rest_chain_triplet_list** (*str*) - (None) Triplet list composed by (reference group, restrain group, chain) list.
29
30 Examples:
31 This is a use example of how to use the building block from Python::
32
33 from biobb_gromacs.gromacs_extra.ndx2resttop import ndx2resttop
34 prop = { 'ref_rest_chain_triplet_list': '( Chain_A, Chain_A_noMut, A ), ( Chain_B, Chain_B_noMut, B ), ( Chain_C, Chain_C_noMut, C ), ( Chain_D, Chain_D_noMut, D )' }
35 ndx2resttop(input_ndx_path='/path/to/myIndex.ndx',
36 input_top_zip_path='/path/to/myTopology.zip',
37 output_top_zip_path='/path/to/newTopology.zip',
38 properties=prop)
39
40 Info:
41 * wrapped_software:
42 * name: In house
43 * license: Apache-2.0
44 * ontology:
45 * name: EDAM
46 * schema: http://edamontology.org/EDAM.owl
47 """
48
49 def __init__(self, input_ndx_path: str, input_top_zip_path: str, output_top_zip_path: str,
50 properties: Optional[dict] = None, **kwargs) -> None:
51 properties = properties or {}
52
53 # Call parent class constructor
54 super().__init__(properties)
55 self.locals_var_dict = locals().copy()
56
57 # Input/Output files
58 self.io_dict = {
59 "in": {"input_ndx_path": input_ndx_path, "input_top_zip_path": input_top_zip_path},
60 "out": {"output_top_zip_path": output_top_zip_path}
61 }
62
63 # Properties specific for BB
64 self.posres_names = properties.get('posres_names')
65 self.force_constants = properties.get('force_constants', '500 500 500')
66 self.ref_rest_chain_triplet_list = properties.get('ref_rest_chain_triplet_list')
67
68 # Check the properties
69 self.check_properties(properties)
70 self.check_arguments()
71
72 @launchlogger
73 def launch(self) -> int:
74 """Execute the :class:`Ndx2resttop <gromacs_extra.ndx2resttop.Ndx2resttop>` object."""
75 # Setup Biobb
76 if self.check_restart():
77 return 0
78
79 top_file = fu.unzip_top(zip_file=self.io_dict['in'].get("input_top_zip_path", ""), out_log=self.out_log)
80
81 # Create group dictionary from the ndx index file
82 # Each group will have a list with its starting and ending lines in the ndx file
83 groups_dic: dict[str, Any] = {}
84 current_group = ''
85 previous_group = ''
- W293 Blank line contains whitespace
87 # Iterate over the lines of the ndx file
88 ndx_lines = open(self.io_dict['in'].get("input_ndx_path", "")).read().splitlines()
89 for index, line in enumerate(ndx_lines):
- W293 Blank line contains whitespace
91 # Found a new group in the index file
92 if line.startswith('['):
93 current_group = line
- W293 Blank line contains whitespace
95 # Set the starting line of the new group
96 groups_dic[current_group] = [index, 0]
97
98 # Check if there was a previous group to set its ending line
99 if previous_group != '':
100 groups_dic[previous_group][1] = index
- W293 Blank line contains whitespace
102 # Update the previous group variable
103 previous_group = current_group
104
105 # If we reach the last line of the file, set the ending line of the last group
106 if index == len(ndx_lines)-1:
107 groups_dic[current_group][1] = index
108
109 # Catch groups with just one line
110 for group_name in groups_dic.keys():
111 if (groups_dic[group_name][0]+1) == groups_dic[group_name][1]:
112 groups_dic[group_name][1] += 1
113
114 fu.log('groups_dic: '+str(groups_dic), self.out_log, self.global_log)
115
116 self.ref_rest_chain_triplet_list = [tuple(elem.strip(' ()').replace(' ', '').split(',')) for elem in str(self.ref_rest_chain_triplet_list).split('),')]
117 fu.log('ref_rest_chain_triplet_list: ' + str(self.ref_rest_chain_triplet_list), self.out_log, self.global_log)
118
119 if self.posres_names:
120 self.posres_names = [elem.strip() for elem in self.posres_names.split()]
121 fu.log('posres_names: ' + str(self.posres_names), self.out_log, self.global_log)
122 else:
123 self.posres_names = ['CUSTOM_POSRES']*len(self.ref_rest_chain_triplet_list)
124 fu.log('posres_names: ' + str(self.posres_names), self.out_log, self.global_log)
125
126 # Check if the number of posres_names matches the number of ref_rest_chain_triplet_list
127 if len(self.posres_names) != len(self.ref_rest_chain_triplet_list):
128 raise ValueError("If posres_names is provided, it should match the number of ref_rest_chain_triplet_list")
129
130 for triplet, posre_name in zip(self.ref_rest_chain_triplet_list, self.posres_names):
131
132 reference_group, restrain_group, chain = triplet
133
134 fu.log('Reference group: '+reference_group, self.out_log, self.global_log)
135 fu.log('Restrain group: '+restrain_group, self.out_log, self.global_log)
136 fu.log('Chain: '+chain, self.out_log, self.global_log)
137 self.io_dict['out']["output_itp_path"] = fu.create_name(path=str(Path(top_file).parent), prefix=self.prefix, step=self.step, name=restrain_group+'_posre.itp')
138
139 # Mapping atoms from absolute enumeration to Chain relative enumeration
140 fu.log('reference_group_index: start_closed:'+str(groups_dic['[ '+reference_group+' ]'][0]+1)+' stop_open: '+str(groups_dic['[ '+reference_group+' ]'][1]), self.out_log, self.global_log)
141 reference_group_list = [int(elem) for line in ndx_lines[groups_dic['[ '+reference_group+' ]'][0]+1: groups_dic['[ '+reference_group+' ]'][1]] for elem in line.split()]
142 fu.log('restrain_group_index: start_closed:'+str(groups_dic['[ '+restrain_group+' ]'][0]+1)+' stop_open: '+str(groups_dic['[ '+restrain_group+' ]'][1]), self.out_log, self.global_log)
143 restrain_group_list = [int(elem) for line in ndx_lines[groups_dic['[ '+restrain_group+' ]'][0]+1: groups_dic['[ '+restrain_group+' ]'][1]] for elem in line.split()]
144 selected_list = [reference_group_list.index(atom)+1 for atom in restrain_group_list]
145 # Creating new ITP with restrains
146 with open(self.io_dict['out'].get("output_itp_path", ''), 'w') as f:
147 fu.log('Creating: '+str(f)+' and adding the selected atoms force constants', self.out_log, self.global_log)
148 f.write('[ position_restraints ]\n')
149 f.write('; atom type fx fy fz\n')
150 for atom in selected_list:
151 f.write(str(atom)+' 1 '+self.force_constants+'\n')
152
153 multi_chain = False
154 # For multi-chain topologies
155 # Including new ITP in the corresponding ITP-chain file
156 for file_dir in Path(top_file).parent.iterdir():
157 if "posre" not in file_dir.name and not file_dir.name.endswith("_pr.itp"):
158 if fnmatch.fnmatch(str(file_dir), "*_chain_"+chain+".itp"):
159 multi_chain = True
160 with open(str(file_dir), 'a') as f:
161 fu.log('Opening: '+str(f)+' and adding the ifdef include statement', self.out_log, self.global_log)
162 f.write('\n')
163 f.write('; Include Position restraint file\n')
164 f.write('#ifdef '+str(posre_name)+'\n')
165 f.write('#include "'+str(Path(self.io_dict['out'].get("output_itp_path", "")).name)+'"\n')
166 f.write('#endif\n')
167 f.write('\n')
168
169 # For single-chain topologies
170 # Including new ITP in the TOP file
171 if not multi_chain:
172
173 # Read all lines of the top file
174 with open(top_file, 'r') as f:
175 lines = f.readlines()
176
177 main_chain = False
178 index = 0
179
180 # Find the index of the line where the custom position restraints are going to be included
181 for line in lines:
182
183 # Find the moleculetype directive of the main chain
184 if line.startswith('Protein_chain_'+chain):
185 main_chain = True
186 index = lines.index(line) + 3
187
188 # Find the end of the moleculetype directive of the main chain
189 if main_chain:
190 if line.startswith('[system]') or line.startswith('[molecules]') or line.startswith('#include ') or line.startswith('#ifdef POSRES'):
191 index = lines.index(line) - 1
192 break
193
194 if index == 0:
195 raise ValueError(f"Protein_chain_{chain} not found in the topology file")
196
197 # Include the custom position restraints in the top file
198 lines.insert(index, '\n')
199 lines.insert(index + 1, '; Include Position restraint file\n')
200 lines.insert(index + 2, '#ifdef '+str(posre_name)+'\n')
201 lines.insert(index + 3, '#include "'+str(Path(self.io_dict['out'].get("output_itp_path", "")).name)+'"\n')
202 lines.insert(index + 4, '#endif\n')
203
204 # Write the new top file
205 with open(top_file, 'w') as f:
206 f.writelines(lines)
207
208 # zip topology
209 fu.zip_top(zip_file=self.io_dict['out'].get("output_top_zip_path", ""), top_file=top_file, out_log=self.out_log, remove_original_files=self.remove_tmp)
210
211 # Remove temporal files
212 self.remove_tmp_files()
213
214 self.check_arguments(output_files_created=True, raise_exception=False)
215 return 0
216
217
218 def ndx2resttop(input_ndx_path: str, input_top_zip_path: str, output_top_zip_path: str,
219 properties: Optional[dict] = None, **kwargs) -> int:
220 """Create :class:`Ndx2resttop <gromacs_extra.ndx2resttop.Ndx2resttop>` class and
221 execute the :meth:`launch() <gromacs_extra.ndx2resttop.Ndx2resttop.launch>` method."""
222 return Ndx2resttop(input_ndx_path=input_ndx_path,
223 input_top_zip_path=input_top_zip_path,
224 output_top_zip_path=output_top_zip_path,
225 properties=properties, **kwargs).launch()
226
227
228 ndx2resttop.__doc__ = Ndx2resttop.__doc__
229
230
231 def main():
232 """Command line execution of this building block. Please check the command line documentation."""
233 parser = argparse.ArgumentParser(description="Wrapper for the GROMACS extra ndx2resttop module.",
234 formatter_class=lambda prog: argparse.RawTextHelpFormatter(prog, width=99999))
235 parser.add_argument('-c', '--config', required=False, help="This file can be a YAML file, JSON file or JSON string")
236
237 # Specific args of each building block
238 required_args = parser.add_argument_group('required arguments')
239 required_args.add_argument('--input_ndx_path', required=True)
240 required_args.add_argument('--input_top_zip_path', required=True)
241 required_args.add_argument('--output_top_zip_path', required=True)
242
243 args = parser.parse_args()
244 config = args.config if args.config else None
245 properties = settings.ConfReader(config=config).get_prop_dic()
246
247 # Specific call of each building block
248 ndx2resttop(input_ndx_path=args.input_ndx_path, input_top_zip_path=args.input_top_zip_path,
249 output_top_zip_path=args.output_top_zip_path, properties=properties)
250
251
252 if __name__ == '__main__':
253 main()