Coverage for mddb_workflow/analyses/rmsf.py: 84%

1# Generic analyses 

2# Easy and fast trajectory analyses carried by Gromacs 

3 

4from numpy import mean, std 

5 

6from mddb_workflow.utils.auxiliar import save_json 

7from mddb_workflow.utils.constants import OUTPUT_RMSF_FILENAME 

8from mddb_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 universe : 'Universe', 

18 output_directory : str, 

19 pbc_selection : 'Selection'): 

20 """Perform the fluctuation analysis.""" 

21 # Set the main output filepath 

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

23 

24 # Run the analysis in all atoms 

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

26 all_atoms = universe.select_atoms('all') 

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

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

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

30 rmsf_values = list(rmsf_analysis.results.rmsf) 

31 

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

33 for index in pbc_selection.atom_indices: 

34 rmsf_values[index] = None 

35 

36 # Get all rmsf values which are not None 

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

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

39 if len(actual_rmsf_values) == 0: 

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

41 return 

42 

43 # Format data 

44 rmsf_data = { 

45 'y': { 

46 'rmsf': { 

47 'average': mean(actual_rmsf_values), 

48 'stddev': std(actual_rmsf_values), 

49 'min': min(actual_rmsf_values), 

50 'max': max(actual_rmsf_values), 

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

52 } 

53 }, 

54 'version': '0.1.0' 

55 } 

56 

57 # Export formatted data to a json file 

58 save_json(rmsf_data, output_analysis_filepath)