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 move_output_file(tmp_out, self.stage_io_dict["out"]["output_ndx_path"], self.out_log, self.global_log) 

134 # Fatslim ignore H atoms so we add them manually 

135 if self.return_hydrogen: 

136 # Parse the atoms indices of the membrane without Hs 

137 leaflet_groups = parse_index(self.stage_io_dict["out"]["output_ndx_path"]) 

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

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

140 # Select the residues using atom indexes 

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

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

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

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

145 # Copy files to host 

146 self.copy_to_host() 

147 

148 # Remove temporary files 

149 self.remove_tmp_files() 

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

151 

152 return self.return_code 

153 

154 

155def 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: 

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

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

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

159 

160 

161fatslim_membranes.__doc__ = FatslimMembranes.__doc__ 

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

163 

164 

165def parse_index(ndx): 

166 """ 

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

168 

169 Args: 

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

171 Returns: 

172 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. 

173 """ 

174 

175 # Read the leaflet.ndx file 

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

177 leaflet_data = file.readlines() 

178 

179 # Initialize dictionaries to store leaflet groups 

180 leaflet_groups = {} 

181 current_group = None 

182 

183 # Parse the leaflet.ndx file 

184 for line in leaflet_data: 

185 line = line.strip() 

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

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

188 leaflet_groups[current_group] = [] 

189 elif current_group is not None: 

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

191 return leaflet_groups 

192 

193 

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

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

196 """ 

197 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. 

198 

199 Args: 

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

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

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

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

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

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

206 Returns: 

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

208 """ 

209 try: 

210 import nglview as nv 

211 except ImportError: 

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

213 

214 u = mda.Universe(topology=input_top_path, 

215 coordinates=input_traj_path) 

216 # Visualize the system with NGLView 

217 view = nv.show_mdanalysis(u) 

218 view.clear_representations() 

219 

220 leaflet_groups = parse_index(output_ndx_path) 

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

222 

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

224 for n in range(n_mems): 

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

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

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

228 non_mem_resn -= set(top_resn) 

229 non_mem_resn -= set(bot_resn) 

230 if leaflets: 

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

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

233 else: 

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

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

236 if len(non_mem_resn) > 0: 

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

238 return view 

239 

240 

241if __name__ == '__main__': 

242 main()