Coverage for model_workflow/analyses/rmsf.py: 86%

1# Generic analyses 

2# Easy and fast trajectory analyses carried by Gromacs 

3 

4from numpy import mean, std 

5 

6from model_workflow.utils.auxiliar import save_json 

7from model_workflow.utils.constants import GREY_HEADER, COLOR_END, OUTPUT_RMSF_FILENAME 

8from model_workflow.utils.type_hints import * 

9 

10from MDAnalysis import Universe 

11from MDAnalysis.analysis import rms 

12 

13 

14# LORE: Fluctuation analysis was done with Gromacs before 

15# LORE: However Gromacs required masses and atom radii, which was a problem in coarse grain 

16def rmsf ( 

17 structure_file : 'File', 

18 trajectory_file : 'File', 

19 output_directory : str, 

20 pbc_selection : 'Selection'): 

21 """Perform the fluctuation analysis.""" 

22 # Set the main output filepath 

23 output_analysis_filepath = f'{output_directory}/{OUTPUT_RMSF_FILENAME}' 

24 

25 # Load the MDAnalysis universe 

26 # Make MDAnalysis warnings and logs grey 

27 print(GREY_HEADER, end='\r') 

28 universe = Universe(structure_file.path, trajectory_file.path) 

29 print(COLOR_END, end='\r') 

30 

31 # Run the analysis in all atoms 

32 # Even atoms to be excluded (e.g. PBC) are excluded later to don't mess atom results order 

33 all_atoms = universe.select_atoms('all') 

34 rmsf_analysis = rms.RMSF(all_atoms).run() 

35 # Make a list of it because it is a ndarray, which is not JSON serializable 

36 # Access to the values thorugh 'results.rmsf' instead of a direct 'rmsf' to avoid a deprecation warning 

37 rmsf_values = list(rmsf_analysis.results.rmsf) 

38 

39 # Filter out values from PBC residue atoms since they may have not sense 

40 for index in pbc_selection.atom_indices: 

41 rmsf_values[index] = None 

42 

43 # Get all rmsf values which are not None 

44 actual_rmsf_values = [ v for v in rmsf_values if v != None ] 

45 # There may be none if the whole system is in PBC 

46 if len(actual_rmsf_values) == 0: 

47 print(' No actual values to do RMSF') 

48 return 

49 

50 # Format data 

51 rmsf_data = { 

52 'y': { 

53 'rmsf': { 

54 'average': mean(actual_rmsf_values), 

55 'stddev': std(actual_rmsf_values), 

56 'min': min(actual_rmsf_values), 

57 'max': max(actual_rmsf_values), 

58 'data': rmsf_values # Keep all values here to make the list length match the number of atoms 

59 } 

60 }, 

61 'version': '0.1.0' 

62 } 

63 

64 # Export formatted data to a json file 

65 save_json(rmsf_data, output_analysis_filepath)