Coverage for model_workflow/analyses/area_per_lipid.py: 100%

1from biobb_mem.fatslim.fatslim_apl import fatslim_apl 

2from model_workflow.utils.auxiliar import save_json 

3from model_workflow.utils.constants import OUTPUT_APL_FILENAME 

4from model_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 } 

37 apl_tmp = '.apl.csv' 

38 print(' Running BioBB FATSLiM APL') 

39 with redirect_stdout(None): 

40 fatslim_apl(input_top_path=structure_file.path, 

41 input_traj_path=trajectory_file.path, 

42 output_csv_path=apl_tmp, 

43 properties=prop) 

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

45 os.remove(apl_tmp) 

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

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

48 # Save the data 

49 data = { 'data':{ 

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

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

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

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

54 'median': m, 

55 'std': s, 

56 } 

57 } 

58 save_json(data, output_analysis_filepath) 

59 

60 

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

62 df = pd.read_csv(output_csv_path) 

63 grids = [] 

64 # Create separate plots for each leaflet 

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

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

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

68 

69 # Define common grid for both plots 

70 x_all = df['X coords'] 

71 y_all = df['Y coords'] 

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

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

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

75 df_leaflet = df[df['leaflet'] == leaflet] 

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

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

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

79 grids.append(grid) 

80 return grids, grid_x, grid_y, m, s