Coverage for model_workflow/analyses/lipid_order.py: 97%
90 statements
« prev ^ index » next coverage.py v7.9.2, created at 2025-07-23 10:54 +0000
1from model_workflow.tools.get_reduced_trajectory import calculate_frame_step
2from model_workflow.utils.mda_spells import to_MDAnalysis_topology
3from model_workflow.utils.auxiliar import save_json
4from model_workflow.utils.constants import OUTPUT_LIPID_ORDER_FILENAME
5from model_workflow.utils.type_hints import *
6import numpy as np
7import MDAnalysis
10def lipid_order (
11 trajectory_file : 'File',
12 standard_topology_file : 'File',
13 output_directory : str,
14 membrane_map: dict,
15 lipid_map: dict,
16 snapshots : int,
17 frames_limit: int = 100):
18 """
19 Calculate lipid order parameters for membranes.
20 This function computes the order parameter (S) for lipid acyl chains, defined as:
21 S = 0.5*(3*<cos²θ> - 1)
22 where θ is the angle between the C-H bond and the membrane normal (z-axis).
24 Notes:
25 - For non-cholesterol lipids, analyzes all acyl chains separately.
26 - For cholesterol, uses carbon atoms bonded to hydrogen.
27 - The order parameter is calculated for each carbon atom in the chain.
28 """
29 if membrane_map is None or membrane_map['n_mems'] == 0:
30 print('-> Skipping lipid order analysis')
31 return
33 # Set the main output filepath
34 output_analysis_filepath = f'{output_directory}/{OUTPUT_LIPID_ORDER_FILENAME}'
36 mda_top = to_MDAnalysis_topology(standard_topology_file)
37 u = MDAnalysis.Universe(mda_top, trajectory_file.path)
38 order_parameters_dict = {}
39 frame_step, _ = calculate_frame_step(snapshots, frames_limit)
40 for ref in lipid_map:
41 # Take the first residue of the reference
42 res = u.residues[ref["resindices"][0]]
43 if ref['lipidmaps']['name'] != 'Cholesterol':
44 carbon_groups = get_all_acyl_chains(res)
45 else:
46 carbon_groups = [res.atoms.select_atoms('element C and bonded element H').indices]
48 order_parameters_dict[ref["resname"]] = {}
49 # For every 'tail'
50 for chain_idx, group in enumerate(carbon_groups):
51 atoms = res.universe.atoms[group]
52 C_names = sorted([atom.name for atom in atoms],key=natural_sort_key)
53 # Find all C-H bonds indices
54 ch_pairs = find_CH_bonds(u, ref["resindices"], C_names)
55 # Initialize the order parameters to sum over the trajectory
56 order_parameters = []
57 costheta_sums = {C_name: np.zeros(len(ch_pairs[C_name]['C'])) for C_name in C_names}
58 n = 0
59 # Loop over the trajectory
60 for ts in u.trajectory[0:snapshots:frame_step]:
61 for C_name in C_names:
62 d = u.atoms[ch_pairs[C_name]['C']].positions - u.atoms[ch_pairs[C_name]['H']].positions
63 costheta_sums[C_name] += d[:,2]**2/np.linalg.norm(d, axis=1)**2
64 n += 1
65 order_parameters = 1.5 * np.array([costheta_sums[C_name].mean() for C_name in C_names])/n - 0.5
66 serrors = 1.5 * np.array([costheta_sums[C_name].std() for C_name in C_names])/n
67 order_parameters_dict[ref["resname"]][str(chain_idx)] = {
68 'atoms': C_names,
69 'avg': order_parameters.tolist(),
70 'std': serrors.tolist()}
71 # Save the data
72 data = { 'data': order_parameters_dict}
73 save_json(data, output_analysis_filepath)
76def get_all_acyl_chains(residue: 'MDAnalysis.Residue') -> list:
77 """
78 Finds all groups of connected Carbon atoms within a residue, including cyclic structures.
80 Parameters
81 ----------
82 residue : MDAnalysis.core.groups.Residue
83 The residue to analyze
85 Returns
86 -------
87 list
88 A list of lists, where each inner list contains atom indices of
89 connected Carbon atoms forming a distinct group
90 """
91 def explore_carbon_group(start_atom, visited):
92 """Helper function to explore a connected group of carbons"""
93 to_visit = [start_atom]
94 group = set()
95 while to_visit:
96 current_atom = to_visit.pop(0)
97 if current_atom.index in visited:
98 continue
99 visited.add(current_atom.index)
100 if current_atom.element == 'C':
101 group.add(current_atom.index)
102 # Add all bonded carbon atoms to the visit queue
103 for bond in current_atom.bonds:
104 for bonded_atom in bond.atoms:
105 if (bonded_atom.element == 'C' and
106 bonded_atom.index not in visited):
107 to_visit.append(bonded_atom)
108 return list(group)
110 # Get all Carbon atoms in the residue
111 carbon_atoms = residue.atoms.select_atoms('element C')
112 visited = set()
113 # Remove ester carbons (without hydrogens)
114 for at in carbon_atoms:
115 if 'H' not in at.bonded_atoms.elements:
116 visited.add(at.index)
117 carbon_groups = []
118 # Find all distinct groups of connected carbons
119 for carbon in carbon_atoms:
120 if carbon.index not in visited:
121 # Get all carbons connected to this one
122 connected_carbons = explore_carbon_group(carbon, visited)
123 if len(connected_carbons)>6: # Only add non-empty groups
124 carbon_groups.append(sorted(connected_carbons))
125 return carbon_groups
128def find_CH_bonds(universe, lipid_resindices, atom_names):
129 """
130 Find all carbon-hydrogen bonds in a lipid molecule using topology information.
132 Args:
133 universe: MDAnalysis Universe object with topology information
134 lipid_resname: Residue name of the lipid (default: 'DPPC')
136 Returns:
137 dict: Dict for each carbon, of tuples containing (carbon_index, hydrogen_index) for each C-H bond
138 """
139 # Get all carbons in the lipid
140 carbons = universe.residues[lipid_resindices].atoms.select_atoms('name ' + ' '.join(atom_names))
141 # Get how many carbons with different name that we will average over later
142 ch_pairs = {}
143 for name in atom_names:
144 ch_pairs[name] = { 'C':[], 'H':[] }
146 # Iterate through each carbon
147 for carbon in carbons:
148 # Get all bonds for this carbon
149 for bond in carbon.bonds:
150 # Check if the bond is between C and H
151 atoms = bond.atoms
152 at1_nm = atoms[0].name
153 at2_nm = atoms[1].name
154 if ('C' in at1_nm and 'H' in at2_nm):
155 # Store the indices making sure carbon is first
156 ch_pairs[carbon.name]['C'].append(atoms[0].index)
157 ch_pairs[carbon.name]['H'].append(atoms[1].index)
158 elif 'H' in at1_nm and 'C' in at2_nm:
159 ch_pairs[carbon.name]['C'].append(atoms[1].index)
160 ch_pairs[carbon.name]['H'].append(atoms[0].index)
161 # Convert lists to numpy arrays
162 for name in ch_pairs:
163 ch_pairs[name]['C'] = np.array(ch_pairs[name]['C'])
164 ch_pairs[name]['H'] = np.array(ch_pairs[name]['H'])
165 return ch_pairs
167def natural_sort_key(s):
168 """Sort strings that contain numbers in human order (C1, C2, C3 instead of C1, C11, C2)"""
169 import re
170 return [int(text) if text.isdigit() else text.lower()
171 for text in re.split('([0-9]+)', s)]