Coverage for mddb_workflow/analyses/rgyr.py: 93%

27 statements  

« prev     ^ index     » next       coverage.py v7.11.0, created at 2025-10-29 15:48 +0000

1from os import remove 

2from numpy import mean, std 

3 

4from mddb_workflow.tools.get_reduced_trajectory import get_reduced_trajectory 

5from mddb_workflow.tools.xvg_parse import xvg_parse 

6from mddb_workflow.utils.auxiliar import save_json 

7from mddb_workflow.utils.constants import OUTPUT_RGYR_FILENAME 

8from mddb_workflow.utils.gmx_spells import run_gromacs 

9from mddb_workflow.utils.type_hints import * 

10 

11# Set an auxiliar data filename 

12rgyr_data_filename = '.rgyr_data.xvg' 

13 

14 

15def rgyr ( 

16 structure_file : 'File', 

17 trajectory_file : 'File', 

18 output_directory : str, 

19 snapshots : int, 

20 structure : 'Structure', 

21 pbc_selection : 'Selection', 

22 frames_limit : int = 5000 

23): 

24 """Perform the RMSd analysis. Use the first trajectory frame in .pdb format as a reference.""" 

25 

26 # Set the main output filepath 

27 output_analysis_filepath = f'{output_directory}/{OUTPUT_RGYR_FILENAME}' 

28 

29 # Use a reduced trajectory in case the original trajectory has many frames 

30 reduced_trajectory_filepath, step, frames = get_reduced_trajectory( 

31 structure_file, 

32 trajectory_file, 

33 snapshots, 

34 frames_limit, 

35 ) 

36 

37 # Generate a custom index file to exclude PBC residues from the analysis 

38 not_pbc_selection = structure.invert_selection(pbc_selection) 

39 if not not_pbc_selection: 

40 print(' No selection to run the analysis') 

41 return 

42 selection_name = 'pbc_atoms' 

43 ndx_selection = not_pbc_selection.to_ndx(selection_name) 

44 ndx_filename = '.not_pbc.ndx' 

45 with open(ndx_filename, 'w') as file: 

46 file.write(ndx_selection) 

47 

48 # Run Gromacs 

49 run_gromacs(f'gyrate -s {structure_file.path} -f {reduced_trajectory_filepath} \ 

50 -o {rgyr_data_filename} -n {ndx_filename}', user_input = selection_name) 

51 

52 # Read the output file and parse it 

53 raw_rgyr_data = xvg_parse(rgyr_data_filename, ['times', 'rgyr', 'rgyrx', 'rgyry', 'rgyrz']) 

54 

55 # Format data 

56 rgyr_data = { 

57 'start': 0, 

58 'step': step, 

59 'y': { 

60 'rgyr': { 

61 'average': mean(raw_rgyr_data['rgyr']), 

62 'stddev': std(raw_rgyr_data['rgyr']), 

63 'min': min(raw_rgyr_data['rgyr']), 

64 'max': max(raw_rgyr_data['rgyr']), 

65 'data': raw_rgyr_data['rgyr'] 

66 }, 

67 'rgyrx': { 

68 'average': mean(raw_rgyr_data['rgyrx']), 

69 'stddev': std(raw_rgyr_data['rgyrx']), 

70 'min': min(raw_rgyr_data['rgyrx']), 

71 'max': max(raw_rgyr_data['rgyrx']), 

72 'data': raw_rgyr_data['rgyrx'] 

73 }, 

74 'rgyry': { 

75 'average': mean(raw_rgyr_data['rgyry']), 

76 'stddev': std(raw_rgyr_data['rgyry']), 

77 'min': min(raw_rgyr_data['rgyry']), 

78 'max': max(raw_rgyr_data['rgyry']), 

79 'data': raw_rgyr_data['rgyry'] 

80 }, 

81 'rgyrz': { 

82 'average': mean(raw_rgyr_data['rgyrz']), 

83 'stddev': std(raw_rgyr_data['rgyrz']), 

84 'min': min(raw_rgyr_data['rgyrz']), 

85 'max': max(raw_rgyr_data['rgyrz']), 

86 'data': raw_rgyr_data['rgyrz'] 

87 } 

88 } 

89 } 

90 

91 # Export formatted data to a json file 

92 save_json(rgyr_data, output_analysis_filepath) 

93 

94 # Remove residual files 

95 remove(rgyr_data_filename) 

96 remove(ndx_filename)