Coverage for biobb_mem/fatslim/fatslim_membranes.py: 56%
105 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-05-28 14:30 +0000
1#!/usr/bin/env python3
3"""Module containing the FATSLiM Membranes class and the command line interface."""
4from pathlib import Path, PurePath
5from biobb_common.generic.biobb_object import BiobbObject
6from biobb_common.tools.file_utils import launchlogger
7from biobb_mem.fatslim.common import ignore_no_box, move_output_file
8import MDAnalysis as mda
9import numpy as np
12class FatslimMembranes(BiobbObject):
13 """
14 | biobb_mem FatslimMembranes
15 | Wrapper of the `FATSLiM membranes <https://pythonhosted.org/fatslim/documentation/leaflets.html>`_ module for leaflet and membrane identification.
16 | FATSLiM is designed to provide efficient and robust analysis of physical parameters from MD trajectories, with a focus on processing large trajectory files quickly.
18 Args:
19 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).
20 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).
21 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).
22 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).
23 properties (dic - Python dictionary object containing the tool parameters, not input/output files):
24 * **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>`_.
25 * **cutoff** (*float*) - (2) Cutoff distance (in nm) to be used when leaflet identification is performed.
26 * **begin_frame** (*int*) - (-1) First frame index to be used for analysis.
27 * **end_frame** (*int*) - (-1) Last frame index to be used for analysis.
28 * **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.
29 * **return_hydrogen** (*bool*) - (False) Include hydrogen atoms in the output index file.
30 * **binary_path** (*str*) - ("fatslim") Path to the fatslim executable binary.
31 * **remove_tmp** (*bool*) - (True) [WF property] Remove temporal files.
32 * **restart** (*bool*) - (False) [WF property] Do not execute if output files exist.
33 * **sandbox_path** (*str*) - ("./") [WF property] Parent path to the sandbox directory.
34 * **container_path** (*str*) - (None) Container path definition.
35 * **container_image** (*str*) - ('afandiadib/ambertools:serial') Container image definition.
36 * **container_volume_path** (*str*) - ('/tmp') Container volume path definition.
37 * **container_working_dir** (*str*) - (None) Container working directory definition.
38 * **container_user_id** (*str*) - (None) Container user_id definition.
39 * **container_shell_path** (*str*) - ('/bin/bash') Path to default shell inside the container.
41 Examples:
42 This is a use example of how to use the building block from Python::
44 from biobb_mem.fatslim.fatslim_membranes import fatslim_membranes
45 prop = {
46 'selection': '(resname DPPC and name P8)',
47 'cutoff': 2.2
48 }
49 fatslim_membranes(input_top_path='/path/to/myTopology.tpr',
50 input_traj_path='/path/to/myTrajectory.xtc',
51 output_ndx_path='/path/to/newIndex.ndx',
52 properties=prop)
54 Info:
55 * wrapped_software:
56 * name: FATSLiM
57 * version: 0.2.2
58 * license: GNU
59 * ontology:
60 * name: EDAM
61 * schema: http://edamontology.org/EDAM.owl
63 """
65 def __init__(self, input_top_path, output_ndx_path, input_traj_path=None,
66 input_ndx_path=None, properties=None, **kwargs) -> None:
67 properties = properties or {}
69 # Call parent class constructor
70 super().__init__(properties)
71 self.locals_var_dict = locals().copy()
73 # Input/Output files
74 self.io_dict = {
75 "in": {"input_top_path": input_top_path,
76 "input_traj_path": input_traj_path,
77 "input_ndx_path": input_ndx_path
78 },
79 "out": {"output_ndx_path": output_ndx_path}
80 }
82 # Properties specific for BB
83 self.selection = properties.get('selection', "not protein and element P")
84 self.cutoff = properties.get('cutoff', 2)
85 self.begin_frame = properties.get('begin_frame', -1)
86 self.end_frame = properties.get('end_frame', -1)
87 self.ignore_no_box = properties.get('ignore_no_box', False)
88 self.return_hydrogen = properties.get('return_hydrogen', False)
89 self.binary_path = properties.get('binary_path', 'fatslim')
90 self.properties = properties
92 # Check the properties
93 self.check_properties(properties)
94 self.check_arguments()
96 @launchlogger
97 def launch(self) -> int:
98 """Execute the :class:`FatslimMembranes <fatslim.fatslim_membranes.FatslimMembranes>` object."""
100 # Setup Biobb
101 if self.check_restart():
102 return 0
103 self.stage_files()
105 unique_dir = self.stage_io_dict.get("unique_dir", "")
106 if self.container_path:
107 working_dir = self.container_volume_path if self.container_volume_path else "/data"
108 else:
109 working_dir = unique_dir
111 host_top_path = str(Path(unique_dir).joinpath(PurePath(self.io_dict["in"]["input_top_path"]).name))
112 host_traj_path = None
113 if self.io_dict["in"].get("input_traj_path"):
114 host_traj_path = str(Path(unique_dir).joinpath(PurePath(self.io_dict["in"]["input_traj_path"]).name))
116 # Create index file using MDAnalysis
117 u = mda.Universe(topology=host_top_path, coordinates=host_traj_path)
118 ignore_no_box(u, self.ignore_no_box, self.out_log, self.global_log)
120 # Build the index to select the atoms from the membrane
121 if self.stage_io_dict["in"].get('input_ndx_path', None):
122 tmp_ndx = self.stage_io_dict["in"]["input_ndx_path"]
123 else:
124 tmp_ndx = str(Path(unique_dir).joinpath('_headgroups.ndx'))
125 with mda.selections.gromacs.SelectionWriter(tmp_ndx, mode='w') as ndx:
126 ndx.write(u.select_atoms(self.selection), name='headgroups')
128 if self.stage_io_dict["in"]["input_top_path"].endswith('gro'):
129 cfg = self.stage_io_dict["in"]["input_top_path"]
130 else:
131 # Convert topology .gro and add box dimensions if not available in the topology
132 cfg = str(Path(unique_dir).joinpath('_output.gro'))
133 # Save as GRO file with box information
134 u.atoms.write(cfg)
136 tmp_out = str(Path(unique_dir).joinpath('_output.ndx'))
137 # Build command
138 self.cmd = [
139 "cd", working_dir, ";", self.binary_path, "membranes",
140 "-n", PurePath(tmp_ndx).name,
141 "-c", PurePath(cfg).name,
142 "--output-index", PurePath(tmp_out).name,
143 "--cutoff", str(self.cutoff),
144 "--begin-frame", str(self.begin_frame),
145 "--end-frame", str(self.end_frame)
146 ]
148 # Run Biobb block
149 self.run_biobb()
150 staged_output_ndx = str(Path(unique_dir).joinpath(PurePath(self.stage_io_dict["out"]["output_ndx_path"]).name))
151 move_output_file(tmp_out, staged_output_ndx, self.out_log, self.global_log)
152 # Fatslim ignore H atoms so we add them manually
153 if self.return_hydrogen:
154 # Parse the atoms indices of the membrane without Hs
155 leaflet_groups = parse_index(staged_output_ndx)
156 with mda.selections.gromacs.SelectionWriter(staged_output_ndx, mode='w') as ndx:
157 for key, value in leaflet_groups.items():
158 # Select the residues using atom indexes
159 res_sele = set(u.atoms[np.array(value)-1].residues.resindices)
160 # Use the rexindex to select all the atoms of the residue
161 sele = f"resindex {' '.join(map(str, res_sele))}"
162 ndx.write(u.select_atoms(sele), name=key)
163 # Copy files to host
164 self.copy_to_host()
166 # Remove temporary files
167 self.remove_tmp_files()
168 self.check_arguments(output_files_created=True, raise_exception=False)
170 return self.return_code
173def 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:
174 """Execute the :class:`FatslimMembranes <fatslim.fatslim_membranes.FatslimMembranes>` class and
175 execute the :meth:`launch() <fatslim.fatslim_membranes.FatslimMembranes.launch>` method."""
176 return FatslimMembranes(**dict(locals())).launch()
179fatslim_membranes.__doc__ = FatslimMembranes.__doc__
180main = FatslimMembranes.get_main(fatslim_membranes, "Calculates the density along an axis of a given cpptraj compatible trajectory.")
183def parse_index(ndx):
184 """
185 Parses a GROMACS index file (.ndx) to extract leaflet groups.
187 Args:
188 ndx (str): Path to the GROMACS index file (.ndx).
189 Returns:
190 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.
191 """
193 # Read the leaflet.ndx file
194 with open(ndx, 'r') as file:
195 leaflet_data = file.readlines()
197 # Initialize dictionaries to store leaflet groups
198 leaflet_groups = {}
199 current_group = None
201 # Parse the leaflet.ndx file
202 for line in leaflet_data:
203 line = line.strip()
204 if line.startswith('[') and line.endswith(']'):
205 current_group = line[1:-1].strip()
206 leaflet_groups[current_group] = []
207 elif current_group is not None:
208 leaflet_groups[current_group].extend(map(int, line.split()))
209 return leaflet_groups
212def display_fatslim(input_top_path: str, lipid_sel: str, input_traj_path: str = None, output_ndx_path="leaflets.ndx", leaflets=True,
213 colors=['blue', 'cyan', 'yellow', 'orange', 'purple', 'magenta'], non_mem_color='red'):
214 """
215 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.
217 Args:
218 input_top_path (str): Path to the input topology file.
219 input_traj_path (str, optional): Path to the input trajectory file. Default is None.
220 output_ndx_path (str, optional): Path to the output index file containing leaflet information. Default is "leaflets.ndx".
221 leaflets (bool, optional): If True, visualize individual leaflets. If False, visualize entire membranes. Default is True.
222 colors (list of str, optional): List of colors to use for visualizing the leaflets or membranes. Default is ['blue', 'cyan', 'yellow', 'orange', 'purple', 'magenta'].
223 non_mem_color (str, optional): Color to use for visualizing lipids not in the membrane. Default is 'red'.
224 Returns:
225 nglview.NGLWidget: An NGLView widget displaying the membrane leaflets.
226 """
227 try:
228 import nglview as nv
229 except ImportError:
230 raise ImportError('Please install the nglview package to visualize the membrane/s.')
232 u = mda.Universe(topology=input_top_path,
233 coordinates=input_traj_path)
234 # Visualize the system with NGLView
235 view = nv.show_mdanalysis(u)
236 view.clear_representations()
238 leaflet_groups = parse_index(output_ndx_path)
239 n_mems = len(leaflet_groups.keys())//2
241 non_mem_resn = set(u.select_atoms(lipid_sel).residues.resnums)
242 for n in range(n_mems):
243 # Convert atoms list to resnums (nglview uses cannot use resindex)
244 top_resn = u.atoms[np.array(leaflet_groups[f'membrane_{n+1}_leaflet_1'])-1].residues.resnums
245 bot_resn = u.atoms[np.array(leaflet_groups[f'membrane_{n+1}_leaflet_2'])-1].residues.resnums
246 non_mem_resn -= set(top_resn)
247 non_mem_resn -= set(bot_resn)
248 if leaflets:
249 view.add_point(selection=", ".join(map(str, top_resn)), color=colors[n*2]) # lipids in top leaflet
250 view.add_point(selection=", ".join(map(str, bot_resn)), color=colors[n*2+1]) # lipids in bot leaflet
251 else:
252 mem_resn = np.concatenate((top_resn, bot_resn))
253 view.add_point(selection=", ".join(map(str, mem_resn)), color=colors[n*2]) # lipids in membrane
254 if len(non_mem_resn) > 0:
255 view.add_point(selection=", ".join(map(str, non_mem_resn)), color=non_mem_color) # lipids without membrane
256 return view
259if __name__ == '__main__':
260 main()