Coverage for mddb_workflow/analyses/area_per_lipid.py: 70%

1from biobb_mem.fatslim.fatslim_apl import fatslim_apl 

2from mddb_workflow.utils.auxiliar import save_json, load_json 

3from mddb_workflow.utils.constants import OUTPUT_APL_FILENAME 

4from mddb_workflow.utils.type_hints import * 

5from scipy.interpolate import griddata 

6from contextlib import redirect_stdout 

7import pandas as pd 

8import numpy as np 

9import os 

10 

11 

12def area_per_lipid ( 

13 structure_file : 'File', 

14 trajectory_file : 'File', 

15 output_directory : str, 

16 membrane_map: dict): 

17 """Area per lipid analysis.""" 

18 

19 if membrane_map is None or membrane_map['n_mems'] == 0: 

20 print('-> Skipping area per lipid analysis') 

21 return 

22 

23 # Set the main output filepath 

24 output_analysis_filepath = f'{output_directory}/{OUTPUT_APL_FILENAME}' 

25 

26 head_sel = [] 

27 for n in range(membrane_map['n_mems']): 

28 head_sel.extend(membrane_map['mems'][str(n)]['polar_atoms']['top']) 

29 head_sel.extend(membrane_map['mems'][str(n)]['polar_atoms']['bot']) 

30 head_sel_mda = 'index ' + " ".join(map(str,(head_sel))) 

31 # Run the analysis on the whole membrane 

32 prop = { 

33 'lipid_selection': head_sel_mda, 

34 'ignore_no_box': True, 

35 'disable_logs': True, 

36 'disable_sandbox': True, 

37 } 

38 apl_tmp = '.apl.csv' 

39 print(' Running BioBB FATSLiM APL') 

40 with redirect_stdout(None): 

41 fatslim_apl(input_top_path=structure_file.path, 

42 input_traj_path=trajectory_file.path, 

43 output_csv_path=apl_tmp, 

44 properties=prop) 

45 grids, grid_x, grid_y, m, s = process_apl(apl_tmp) 

46 os.remove(apl_tmp) 

47 # Replace NaNs with -1 in the grids so the loader don't break 

48 grids = [np.nan_to_num(grid, nan=-1) for grid in grids] 

49 # Save the data 

50 data = { 'data':{ 

51 'lower leaflet': grids[0].tolist(), 

52 'upper leaflet': grids[1].tolist(), 

53 'grid_x': grid_x[:,0].tolist(), 

54 'grid_y': grid_y[0,:].tolist(), 

55 'median': m, 

56 'std': s, 

57 } 

58 } 

59 save_json(data, output_analysis_filepath) 

60 

61 

62def process_apl(output_csv_path, res=100j): 

63 df = pd.read_csv(output_csv_path) 

64 grids = [] 

65 # Create separate plots for each leaflet 

66 df['Area per lipid'] *= 100 # Convert to A^2 

67 m = df['Area per lipid'].median() 

68 s = df['Area per lipid'].std() 

69 

70 # Define common grid for both plots 

71 x_all = df['X coords'] 

72 y_all = df['Y coords'] 

73 grid_x, grid_y = np.mgrid[min(x_all):max(x_all):res, 

74 min(y_all):max(y_all):res] 

75 for leaflet in ['lower leaflet', 'upper leaflet']: 

76 df_leaflet = df[df['leaflet'] == leaflet] 

77 points = np.stack((np.array(df_leaflet['X coords']).T, np.array(df_leaflet['Y coords']).T), axis=-1) 

78 values = np.array(df_leaflet['Area per lipid']) 

79 grid = griddata(points, values, (grid_x, grid_y), method='cubic') 

80 grids.append(grid) 

81 return grids, grid_x, grid_y, m, s 

82 

83 

84def plot_apl(output_analysis_filepath): 

85 """Plot the area per lipid for both leaflets.""" 

86 import matplotlib.pyplot as plt 

87 data = load_json(output_analysis_filepath) 

88 grids = data['data'] 

89 lower_leaflet = np.array(grids['lower leaflet']) 

90 upper_leaflet = np.array(grids['upper leaflet']) 

91 grid_x = np.array(grids['grid_x']) 

92 grid_y = np.array(grids['grid_y']) 

93 

94 fig, axs = plt.subplots(1, 2, figsize=(12, 6)) 

95 c1 = axs[0].contourf(grid_x, grid_y, lower_leaflet.T, cmap='viridis') 

96 axs[0].set_title('Lower Leaflet') 

97 axs[0].set_xlabel('X (Å)') 

98 axs[0].set_ylabel('Y (Å)') 

99 fig.colorbar(c1, ax=axs[0]) 

100 

101 c2 = axs[1].contourf(grid_x, grid_y, upper_leaflet.T, cmap='viridis') 

102 axs[1].set_title('Upper Leaflet') 

103 axs[1].set_xlabel('X (Å)') 

104 axs[1].set_ylabel('Y (Å)') 

105 fig.colorbar(c2, ax=axs[1]) 

106 

107 plt.tight_layout() 

108 plt.show()