Coverage for model_workflow/analyses/rgyr.py: 83%

1from os.path import exists 

2from os import remove 

3from subprocess import run, PIPE, Popen 

4from numpy import mean, std 

5 

6from model_workflow.tools.get_reduced_trajectory import get_reduced_trajectory 

7from model_workflow.tools.xvg_parse import xvg_parse 

8from model_workflow.utils.auxiliar import save_json 

9from model_workflow.utils.constants import GROMACS_EXECUTABLE, OUTPUT_RGYR_FILENAME 

10from model_workflow.utils.type_hints import * 

11 

12# Set an auxiliar data filename 

13rgyr_data_filename = '.rgyr_data.xvg' 

14 

15 

16def rgyr ( 

17 structure_file : 'File', 

18 trajectory_file : 'File', 

19 output_directory : str, 

20 snapshots : int, 

21 structure : 'Structure', 

22 pbc_selection : 'Selection', 

23 frames_limit : int = 5000 

24): 

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

26 

27 # Set the main output filepath 

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

29 

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

31 reduced_trajectory_filepath, step, frames = get_reduced_trajectory( 

32 structure_file, 

33 trajectory_file, 

34 snapshots, 

35 frames_limit, 

36 ) 

37 

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

39 not_pbc_selection = structure.invert_selection(pbc_selection) 

40 if not not_pbc_selection: 

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

42 return 

43 selection_name = 'pbc_atoms' 

44 ndx_selection = not_pbc_selection.to_ndx(selection_name) 

45 ndx_filename = '.not_pbc.ndx' 

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

47 file.write(ndx_selection) 

48 

49 # Run Gromacs 

50 p = Popen([ 

51 "echo", 

52 selection_name, 

53 ], stdout=PIPE) 

54 process = run([ 

55 GROMACS_EXECUTABLE, 

56 "gyrate", 

57 "-s", 

58 structure_file.path, 

59 "-f", 

60 reduced_trajectory_filepath, 

61 '-o', 

62 rgyr_data_filename, 

63 '-n', 

64 ndx_filename, 

65 '-quiet' 

66 ], stdin=p.stdout, stdout=PIPE, stderr=PIPE) 

67 logs = process.stdout.decode() 

68 p.stdout.close() 

69 

70 # If the output does not exist at this point it means something went wrong with gromacs 

71 if not exists(rgyr_data_filename): 

72 print(logs) 

73 error_logs = process.stderr.decode() 

74 print(error_logs) 

75 raise SystemExit('Something went wrong with GROMACS') 

76 

77 # Read the output file and parse it 

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

79 

80 # Format data 

81 rgyr_data = { 

82 'start': 0, 

83 'step': step, 

84 'y': { 

85 'rgyr': { 

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

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

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

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

90 'data': raw_rgyr_data['rgyr'] 

91 }, 

92 'rgyrx': { 

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

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

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

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

97 'data': raw_rgyr_data['rgyrx'] 

98 }, 

99 'rgyry': { 

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

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

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

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

104 'data': raw_rgyr_data['rgyry'] 

105 }, 

106 'rgyrz': { 

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

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

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

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

111 'data': raw_rgyr_data['rgyrz'] 

112 } 

113 } 

114 } 

115 

116 # Export formatted data to a json file 

117 save_json(rgyr_data, output_analysis_filepath) 

118 

119 # Remove residual files 

120 remove(rgyr_data_filename) 

121 remove(ndx_filename)