Coverage for biobb_mem/fatslim/fatslim_membranes.py: 53%

1#!/usr/bin/env python3 

2 

3"""Module containing the FATSLiM Membranes class and the command line interface.""" 

4from biobb_common.generic.biobb_object import BiobbObject 

5from biobb_common.tools.file_utils import launchlogger 

6from biobb_mem.fatslim.common import ignore_no_box, move_output_file 

7import MDAnalysis as mda 

8import numpy as np 

9 

10 

11class FatslimMembranes(BiobbObject): 

12 """ 

13 | biobb_mem FatslimMembranes 

14 | Wrapper of the `FATSLiM membranes <https://pythonhosted.org/fatslim/documentation/leaflets.html>`_ module for leaflet and membrane identification. 

15 | FATSLiM is designed to provide efficient and robust analysis of physical parameters from MD trajectories, with a focus on processing large trajectory files quickly. 

16 

17 Args: 

18 input_top_path (str): Path to the input topology file. File type: input. `Sample file <https://github.com/bioexcel/biobb_mem/raw/main/biobb_mem/test/data/A01JD/A01JD.pdb>`_. Accepted formats: tpr (edam:format_2333), gro (edam:format_2033), g96 (edam:format_2033), pdb (edam:format_1476), brk (edam:format_2033), ent (edam:format_1476). 

19 input_traj_path (str) (Optional): Path to the GROMACS trajectory file. File type: input. `Sample file <https://github.com/bioexcel/biobb_mem/raw/main/biobb_mem/test/data/A01JD/A01JD.xtc>`_. Accepted formats: xtc (edam:format_3875), trr (edam:format_3910), cpt (edam:format_2333), gro (edam:format_2033), g96 (edam:format_2033), pdb (edam:format_1476), tng (edam:format_3876). 

20 input_ndx_path (str) (Optional): Path to the input lipid headgroups index NDX file. File type: input. `Sample file <https://github.com/bioexcel/biobb_mem/raw/main/biobb_mem/test/data/A01JD/A01JD.ndx>`_. Accepted formats: ndx (edam:format_2033). 

21 output_ndx_path (str): Path to the output index NDX file. File type: output. `Sample file <https://github.com/bioexcel/biobb_mem/raw/main/biobb_mem/test/reference/fatslim/leaflets.ndx>`_. Accepted formats: ndx (edam:format_2033). 

22 properties (dic - Python dictionary object containing the tool parameters, not input/output files): 

23 * **selection** (*str*) - ("not protein and element P") Alternative ot the NDX file for choosing the Headgroups used in the identification using MDAnalysis `selection language <https://docs.mdanalysis.org/stable/documentation_pages/selections.html>`_. 

24 * **cutoff** (*float*) - (2) Cutoff distance (in nm) to be used when leaflet identification is performed. 

25 * **begin_frame** (*int*) - (-1) First frame index to be used for analysis. 

26 * **end_frame** (*int*) - (-1) Last frame index to be used for analysis. 

27 * **ignore_no_box** (*bool*) - (False) Ignore the absence of box information in the topology. If the topology does not contain box information, the box will be set to the minimum and maximum positions of the atoms. 

28 * **return_hydrogen** (*bool*) - (False) Include hydrogen atoms in the output index file. 

29 * **binary_path** (*str*) - ("fatslim") Path to the fatslim executable binary. 

30 * **remove_tmp** (*bool*) - (True) [WF property] Remove temporal files. 

31 * **restart** (*bool*) - (False) [WF property] Do not execute if output files exist. 

32 * **sandbox_path** (*str*) - ("./") [WF property] Parent path to the sandbox directory. 

33 

34 Examples: 

35 This is a use example of how to use the building block from Python:: 

36 

37 from biobb_mem.fatslim.fatslim_membranes import fatslim_membranes 

38 prop = { 

39 'selection': '(resname DPPC and name P8)', 

40 'cutoff': 2.2 

41 } 

42 fatslim_membranes(input_top_path='/path/to/myTopology.tpr', 

43 input_traj_path='/path/to/myTrajectory.xtc', 

44 output_ndx_path='/path/to/newIndex.ndx', 

45 properties=prop) 

46 

47 Info: 

48 * wrapped_software: 

49 * name: FATSLiM 

50 * version: 0.2.2 

51 * license: GNU 

52 * ontology: 

53 * name: EDAM 

54 * schema: http://edamontology.org/EDAM.owl 

55 

56 """ 

57 

58 def __init__(self, input_top_path, output_ndx_path, input_traj_path=None, 

59 input_ndx_path=None, properties=None, **kwargs) -> None: 

60 properties = properties or {} 

61 

62 # Call parent class constructor 

63 super().__init__(properties) 

64 self.locals_var_dict = locals().copy() 

65 

66 # Input/Output files 

67 self.io_dict = { 

68 "in": {"input_top_path": input_top_path, 

69 "input_traj_path": input_traj_path, 

70 "input_ndx_path": input_ndx_path 

71 }, 

72 "out": {"output_ndx_path": output_ndx_path} 

73 } 

74 

75 # Properties specific for BB 

76 self.selection = properties.get('selection', "not protein and element P") 

77 self.cutoff = properties.get('cutoff', 2) 

78 self.begin_frame = properties.get('begin_frame', -1) 

79 self.end_frame = properties.get('end_frame', -1) 

80 self.ignore_no_box = properties.get('ignore_no_box', False) 

81 self.return_hydrogen = properties.get('return_hydrogen', False) 

82 self.binary_path = properties.get('binary_path', 'fatslim') 

83 self.properties = properties 

84 

85 # Check the properties 

86 self.check_properties(properties) 

87 self.check_arguments() 

88 

89 @launchlogger 

90 def launch(self) -> int: 

91 """Execute the :class:`FatslimMembranes <fatslim.fatslim_membranes.FatslimMembranes>` object.""" 

92 

93 # Setup Biobb 

94 if self.check_restart(): 

95 return 0 

96 self.stage_files() 

97 

98 # Create index file using MDAnalysis 

99 u = mda.Universe(topology=self.stage_io_dict["in"]["input_top_path"], 

100 coordinates=self.stage_io_dict["in"].get("input_traj_path")) 

101 ignore_no_box(u, self.ignore_no_box, self.out_log, self.global_log) 

102 

103 # Build the index to select the atoms from the membrane 

104 if self.stage_io_dict["in"].get('input_ndx_path', None): 

105 tmp_ndx = self.stage_io_dict["in"]["input_ndx_path"] 

106 else: 

107 tmp_ndx = self.create_tmp_file('_headgroups.ndx') 

108 with mda.selections.gromacs.SelectionWriter(tmp_ndx, mode='w') as ndx: 

109 ndx.write(u.select_atoms(self.selection), name='headgroups') 

110 

111 if self.stage_io_dict["in"]["input_top_path"].endswith('gro'): 

112 cfg = self.stage_io_dict["in"]["input_top_path"] 

113 else: 

114 # Convert topology .gro and add box dimensions if not available in the topology 

115 cfg = self.create_tmp_file('_output.gro') 

116 # Save as GRO file with box information 

117 u.atoms.write(cfg) 

118 

119 tmp_out = self.create_tmp_file('_output.ndx') 

120 # Build command 

121 self.cmd = [ 

122 self.binary_path, "membranes", 

123 "-n", tmp_ndx, 

124 "-c", cfg, 

125 "--output-index", tmp_out, 

126 "--cutoff", str(self.cutoff), 

127 "--begin-frame", str(self.begin_frame), 

128 "--end-frame", str(self.end_frame) 

129 ] 

130 

131 # Run Biobb block 

132 self.run_biobb() 

133 # Fatslim ignore H atoms so we add them manually 

134 if self.return_hydrogen: 

135 # Parse the atoms indices of the membrane without Hs 

136 leaflet_groups = parse_index(tmp_out[:-4]+'_0000.ndx') 

137 with mda.selections.gromacs.SelectionWriter(self.stage_io_dict["out"]["output_ndx_path"], mode='w') as ndx: 

138 for key, value in leaflet_groups.items(): 

139 # Select the residues using atom indexes 

140 res_sele = set(u.atoms[np.array(value)-1].residues.resindices) 

141 # Use the rexindex to select all the atoms of the residue 

142 sele = f"resindex {' '.join(map(str, res_sele))}" 

143 ndx.write(u.select_atoms(sele), name=key) 

144 else: 

145 move_output_file(tmp_out, self.stage_io_dict["out"]["output_ndx_path"], self.out_log, self.global_log) 

146 # Copy files to host 

147 self.copy_to_host() 

148 

149 # Remove temporary files 

150 self.remove_tmp_files() 

151 self.check_arguments(output_files_created=True, raise_exception=False) 

152 

153 return self.return_code 

154 

155 

156def fatslim_membranes(input_top_path: str, output_ndx_path: str, input_traj_path: str = None, input_ndx_path: str = None, properties: dict = None, **kwargs) -> int: 

157 """Execute the :class:`FatslimMembranes <fatslim.fatslim_membranes.FatslimMembranes>` class and 

158 execute the :meth:`launch() <fatslim.fatslim_membranes.FatslimMembranes.launch>` method.""" 

159 return FatslimMembranes(**dict(locals())).launch() 

160 

161 

162fatslim_membranes.__doc__ = FatslimMembranes.__doc__ 

163main = FatslimMembranes.get_main(fatslim_membranes, "Calculates the density along an axis of a given cpptraj compatible trajectory.") 

164 

165 

166def parse_index(ndx): 

167 """ 

168 Parses a GROMACS index file (.ndx) to extract leaflet groups. 

169 

170 Args: 

171 ndx (str): Path to the GROMACS index file (.ndx). 

172 Returns: 

173 dict: A dictionary where keys are group names for each leaflet in format "membrane_1_leaflet_1" and values are lists of integers representing atom indices starting from 1. 

174 """ 

175 

176 # Read the leaflet.ndx file 

177 with open(ndx, 'r') as file: 

178 leaflet_data = file.readlines() 

179 

180 # Initialize dictionaries to store leaflet groups 

181 leaflet_groups = {} 

182 current_group = None 

183 

184 # Parse the leaflet.ndx file 

185 for line in leaflet_data: 

186 line = line.strip() 

187 if line.startswith('[') and line.endswith(']'): 

188 current_group = line[1:-1].strip() 

189 leaflet_groups[current_group] = [] 

190 elif current_group is not None: 

191 leaflet_groups[current_group].extend(map(int, line.split())) 

192 return leaflet_groups 

193 

194 

195def display_fatslim(input_top_path: str, lipid_sel: str, input_traj_path: str = None, output_ndx_path="leaflets.ndx", leaflets=True, 

196 colors=['blue', 'cyan', 'yellow', 'orange', 'purple', 'magenta'], non_mem_color='red'): 

197 """ 

198 Visualize the leaflets of a membrane using NGLView. The lipids in the membrane are colored according to their leaflet. The ones not in the membrane are colored in red. 

199 

200 Args: 

201 input_top_path (str): Path to the input topology file. 

202 input_traj_path (str, optional): Path to the input trajectory file. Default is None. 

203 output_ndx_path (str, optional): Path to the output index file containing leaflet information. Default is "leaflets.ndx". 

204 leaflets (bool, optional): If True, visualize individual leaflets. If False, visualize entire membranes. Default is True. 

205 colors (list of str, optional): List of colors to use for visualizing the leaflets or membranes. Default is ['blue', 'cyan', 'yellow', 'orange', 'purple', 'magenta']. 

206 non_mem_color (str, optional): Color to use for visualizing lipids not in the membrane. Default is 'red'. 

207 Returns: 

208 nglview.NGLWidget: An NGLView widget displaying the membrane leaflets. 

209 """ 

210 try: 

211 import nglview as nv 

212 except ImportError: 

213 raise ImportError('Please install the nglview package to visualize the membrane/s.') 

214 

215 u = mda.Universe(topology=input_top_path, 

216 coordinates=input_traj_path) 

217 # Visualize the system with NGLView 

218 view = nv.show_mdanalysis(u) 

219 view.clear_representations() 

220 

221 leaflet_groups = parse_index(output_ndx_path) 

222 n_mems = len(leaflet_groups.keys())//2 

223 

224 non_mem_resn = set(u.select_atoms(lipid_sel).residues.resnums) 

225 for n in range(n_mems): 

226 # Convert atoms list to resnums (nglview uses cannot use resindex) 

227 top_resn = u.atoms[np.array(leaflet_groups[f'membrane_{n+1}_leaflet_1'])-1].residues.resnums 

228 bot_resn = u.atoms[np.array(leaflet_groups[f'membrane_{n+1}_leaflet_2'])-1].residues.resnums 

229 non_mem_resn -= set(top_resn) 

230 non_mem_resn -= set(bot_resn) 

231 if leaflets: 

232 view.add_point(selection=", ".join(map(str, top_resn)), color=colors[n*2]) # lipids in top leaflet 

233 view.add_point(selection=", ".join(map(str, bot_resn)), color=colors[n*2+1]) # lipids in bot leaflet 

234 else: 

235 mem_resn = np.concatenate((top_resn, bot_resn)) 

236 view.add_point(selection=", ".join(map(str, mem_resn)), color=colors[n*2]) # lipids in membrane 

237 if len(non_mem_resn) > 0: 

238 view.add_point(selection=", ".join(map(str, non_mem_resn)), color=non_mem_color) # lipids without membrane 

239 return view 

240 

241 

242if __name__ == '__main__': 

243 main()