Water bridges detection in a trajectory

Now, we will perform calculations for a trajectory file that was obtained using the NAMD_ package. We will use calcWaterBridgesTrajectory(), for which we need to provide PDB and DCD files.

The system (protein in a water box) can be found in 5kqm_all_sci.pdb. The trajectory, NAMD_D2_sample.dcd, has dcd format. If we want to analyze trajectories with different formats, we need to convert them to dcd file format or save the trajectory as a multi-model PDB (using VMD or another tool).

Parse structure with trajectory

In [1]: PDBtraj_file = "5kqm_all_sci.pdb"

In [2]: coords_traj = parsePDB(PDBtraj_file)

In [3]: trajectory = parseDCD("NAMD_D2_sample.dcd")
@> 19321 atoms and 1 coordinate set(s) were parsed in 0.18s.
@> DCD file contains 17 coordinate sets for 19321 atoms.
@> DCD file was parsed in 0.01 seconds.
@> 3.76 MB parsed at input rate 748.06 MB/s.
@> 17 coordinate sets parsed at input rate 3382 frame/s.

The analysis od water bridges can be performed on selected frames by using start_frame or stop_frame.

In [4]: wb_traj = calcWaterBridgesTrajectory(coords_traj, trajectory, start_frame=5,
   ...:                                      stop_frame=15, output='info')
   ...: 
@> Frame: 5
@> 101 water bridges detected.
@> Frame: 6
@> 107 water bridges detected.
@> Frame: 7
@> 90 water bridges detected.
@> Frame: 8
@> 97 water bridges detected.
@> Frame: 9
@> 122 water bridges detected.
@> Frame: 10
@> 101 water bridges detected.
@> Frame: 11
@> 130 water bridges detected.
@> Frame: 12
@> 132 water bridges detected.
@> Frame: 13
@> 126 water bridges detected.
@> Frame: 14
@> 88 water bridges detected.
@> Frame: 15
@> 105 water bridges detected.

Because of the amount of data, detailed results will not be displayed. We instead access the raw data by using output='info'.

In [5]: wb_traj
[[['THR5',
   'OG1_8',
   'P',
   'TRP39',
   'NE1_547',
   'P',
   3.261452627054999,
   1,
   ['3W_13313']],
  ['THR5',
   'O_15',
   'P',
   'ASP86',
   'OD1_1269',
   'P',
   5.986350034086454,
   2,
   ['3W_12974', '3W_18431']],
  ['THR5',
   'O_15',
   'P',
   'LYS110',
   'NZ_1667',
   'P',
   7.375256709599827,
   2,
   ['3W_12974', '3W_18431']],
  ['THR5',
   'O_15',
   'P',
   'LYS6',
   'NZ_32',
   'P',
   6.414308925017051,
   2,
   ['3W_12974', '3W_12152']],
  ['LYS6',
   'NZ_32',
   'P',
   'TYR87',
   'OH_1286',
   'P',
   4.891713264838611,
   1,
   ['3W_9209']]
   ...
   ...
   ]]

Save the results

The results can be saved using saveWaterBridges() in two formats. The txt file will contain all the results for analysis and can be visualized in a text editor, and the wb file will restore data for further analysis. It can be loaded using parseWaterBridges() as shown below.

First, we have to return the calculation without output='info'.

We can suppress the logged output using confProDy() to set the verbosity to 'none'.

In [6]: confProDy(verbosity='none')

In [7]: wb_traj = calcWaterBridgesTrajectory(coords_traj, trajectory,
   ...:                                      stop_frame=15)
   ...: 
In [8]: saveWaterBridges(wb_traj,'wb_saved.txt')

In [9]: saveWaterBridges(wb_traj,'wb_saved.wb')

To load the wb file, use parseWaterBridges() and protein coordinates as follows:

In [10]: waterBridges = parseWaterBridges('wb_saved.wb', coords_traj)

Loaded results from a .wb file are Atomic type and therefore can be used for analysis later.

Analysis of the results:

Information about residues contributing to water bridges

The data can be analyzed using calcWaterBridgesStatistics(). The following analysis provides details about the pairs of residues engaged in water bridges, their frequency of occurrence, and the average distance between them. The standard deviation offers insights into the variation in distance throughout the simulation. Moreover, the analysis can be saved using the filename option.

We can recover logged output using confProDy() again with a different verbosity.

In [11]: confProDy(verbosity='debug')

In [12]: analysisAtomic = calcWaterBridgesStatistics(waterBridges, trajectory,
   ....:                                             filename='data.txt')
   ....: 
@> RES1           RES2           PERC      DIST_AVG  DIST_STD
@> ARG40P         SER7P          12.500    4.901     0.000
@> ASP92P         ARG18P         68.750    4.285     1.159
@> ASN95P         ARG18P         68.750    5.099     1.192
@> GLU23P         PRO20P         12.500    4.571     0.000
@> HSE72P         GLU23P         12.500    3.669     0.458
@> VAL41P         ARG27P         56.250    5.565     0.781
@> SER71P         ARG27P         75.000    6.116     0.445
@> ASN34P         ASP32P         25.000    4.218     0.652
@> GLU37P         SER36P         75.000    3.700     1.154
@> THR84P         ARG40P         50.000    4.235     0.671
@> ARG75P         ASP42P         68.750    3.159     0.652
@> ASN95P         THR46P         62.500    4.067     0.842
@> TYR49P         SER47P         50.000    4.320     0.757
..
..

The output is a dictionary, so we can use dict.items() to inspect it.

In [13]: for item in list(analysisAtomic.items())[:5]:
   ....:    print(item)
   ....: 
((40, 7), {'percentage': 12.5, 'distAvg': 4.9006157, 'distStd': 0.0})
((7, 40), {'percentage': 12.5, 'distAvg': 4.9006157, 'distStd': 0.0})
((92, 18), {'percentage': 68.75, 'distAvg': 4.2853837, 'distStd': 1.159262})
((18, 92), {'percentage': 68.75, 'distAvg': 4.2853837, 'distStd': 1.159262})
((95, 18), {'percentage': 68.75, 'distAvg': 5.0986476, 'distStd': 1.1916962})

To have easier access to the data, we can use getWaterBridgeStatInfo().

In [14]: wb_stat_info = getWaterBridgeStatInfo(analysisAtomic, coords_traj)

In [15]: wb_stat_info
{('SER7P', 'ARG40P'): {'percentage': 12.5,
  'distAvg': 4.9006157,
  'distStd': 0.0},
 ('ARG18P', 'ASP92P'): {'percentage': 68.75,
  'distAvg': 4.2853837,
  'distStd': 1.159262},
 ('ARG18P', 'ASN95P'): {'percentage': 68.75,
  'distAvg': 5.0986476,
  'distStd': 1.1916962},
 ('PRO20P', 'GLU23P'): {'percentage': 12.5,
  'distAvg': 4.571081,
  'distStd': 0.0},
   ...
   ...

To obtain maps of interactions for the protein structure, we can use showWaterBridgeMatrix(), which is equipped with three paramaters: 'percentage' (how often two residues were forming water bridges), 'distAvg' (how close there were on average), and 'distStd' (how stable that water bridge was).

In [16]: showWaterBridgeMatrix(analysisAtomic, 'percentage')
../../_images/traj_percentage.png
In [17]: showWaterBridgeMatrix(analysisAtomic, 'distAvg')
../../_images/traj_distAvg.png
In [18]: showWaterBridgeMatrix(analysisAtomic, 'distStd')
../../_images/traj_distStd.png

Raw data of the matrices can be obtained with calcWaterBridgeMatrix(). The type of the data in the matrix can be selected using the following strings for the second argument: 'percentage', 'distAvg', 'distStd'.

In [19]: M1 = calcWaterBridgeMatrix(analysisAtomic, 'percentage')

In [20]: M2 = calcWaterBridgeMatrix(analysisAtomic, 'distAvg')

In [21]: M3 = calcWaterBridgeMatrix(analysisAtomic, 'distStd')
In [22]: M1
array([[ 0.  ,  0.  ,  0.  , ...,  0.  ,  0.  ,  0.  ],
       [ 0.  ,  0.  ,  0.  , ...,  0.  ,  0.  ,  0.  ],
       [ 0.  ,  0.  ,  0.  , ...,  0.  ,  0.  ,  0.  ],
       ...,
       [ 0.  ,  0.  ,  0.  , ...,  0.  , 12.5 , 31.25],
       [ 0.  ,  0.  ,  0.  , ..., 12.5 ,  0.  , 12.5 ],
       [ 0.  ,  0.  ,  0.  , ..., 31.25, 12.5 ,  0.  ]])
In [23]: M2
array([[0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       ...,
       [0.        , 0.        , 0.        , ..., 0.        , 4.58851337,
        5.82083416],
       [0.        , 0.        , 0.        , ..., 4.58851337, 0.        ,
        3.52366138],
       [0.        , 0.        , 0.        , ..., 5.82083416, 3.52366138,
        0.        ]])
In [24]: M3
array([[0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       ...,
       [0.        , 0.        , 0.        , ..., 0.        , 1.71697354,
        1.38650537],
       [0.        , 0.        , 0.        , ..., 1.71697354, 0.        ,
        1.27207112],
       [0.        , 0.        , 0.        , ..., 1.38650537, 1.27207112,
        0.        ]])

Statistical analysis for water bridges

To visualize the results in a more accessible way, we can use the calcBridgingResiduesHistogram() function, which will show how often each residue was contributing to the water bridges in the trajectory.

In [25]: wb_res_hist = calcBridgingResiduesHistogram(waterBridges)

In [26]: wb_res_hist
../../_images/traj_res_hist.png
[('LEU96P', 1),
 ('MET63P', 1),
 ('PHE152P', 1),
 ('LEU29P', 1),
 ('PRO130P', 1),
 ('PHE85P', 1),
 ('PRO54P', 1),
 ('ILE16P', 1),
 ('CYS148P', 1),
 ('VAL25P', 1),
 ('ILE77P', 1),
 .
 .
 ('ARG75P', 15),
 ('ARG18P', 15),
 ('ARG65P', 15),
 ('ARG40P', 16),
 ('ARG147P', 16),
 ('ARG58P', 16),
 ('ARG27P', 16),
 ('ASP92P', 16),
 ('TYR49P', 16),
 ('LYS102P', 16),
 ('ARG150P', 16),
 ('SER36P', 16)]

The clip option can be used to include different number of results on the histogram.

In [27]: calcBridgingResiduesHistogram(waterBridges, clip=25)
../../_images/traj_res_hist2.png

If we are interested in one particular residue, we can also use calcWaterBridgesDistribution() to find their partners in water bridges. Below we can see results for arginine 147 or aspartic acid 92 from chain P using the nomenclature for them corresponding to the keys of the dictionary.

In [28]: calcWaterBridgesDistribution(waterBridges, 'ARG147P')
[('GLN122P', 8),
 ('ARG150P', 7),
 ('GLN143P', 6),
 ('LYS123P', 6),
 ('GLN124P', 5),
 ('ASP120P', 5),
 ('GLN144P', 3),
 ('THR140P', 2)]
In [29]: calcWaterBridgesDistribution(waterBridges, 'ASP92P')
[('ARG18P', 11),
 ('ASN95P', 10),
 ('SER94P', 5),
 ('MET91P', 5),
 ('ASP129P', 4),
 ('LEU13P', 3),
 ('CYS90P', 1)]

Once we select a pair of residues which are supported by interactions with water molecules, we can use calcWaterBridgesDistribution() to obtain histograms with results such as distances between them (metric='distance'), the number of water molecules which were involved (metric='waters'), and information about residue part which was involved in water bridges, i.e. backbone or side chain (metric='location').

In [30]: calcWaterBridgesDistribution(waterBridges,  'ASP92P', 'ARG18P',
   ....:                              trajectory=trajectory, metric='distance')
   ....: 
../../_images/traj_distribution.png
[5.3736005,
 5.3736005,
 5.167575,
 2.681302,
 5.371548,
 2.6318514,
 3.0394073,
 4.0884595,
 5.4406505,
 3.4112484,
 2.805657,
 5.4176636,
 3.5104342,
 5.991175,
 5.470093,
 3.4345005,
 3.6427624]
In [31]: calcWaterBridgesDistribution(waterBridges, 'ARG147P', 'GLN122P',
   ....:                                                      metric='waters')
   ....: 
../../_images/traj_distribution2.png
[2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2]
In [32]: calcWaterBridgesDistribution(waterBridges, 'ARG147P', 'GLN122P',
   ....:                              trajectory=trajectory, metric='location')
   ....: 
{'ARG147P': {'backbone': 7, 'side': 86},
'GLN122P': {'backbone': 21, 'side': 25}}

Save results as PDB file

The results can be stored as a PDB file using savePDBWaterBridges() (single PDB file, single frame) or using savePDBWaterBridgesTrajectory() to save all the results (large number of frames saved each independently).

5kqm_all_sci_multi_0.pdb 5kqm_all_sci_multi_1.pdb .. 5kqm_all_sci_multi_15.pdb

Those results can be displayed in any program for visualization. The results for the protein structure will be storage in the B-factor (*beta*) column (average values of contributions of each residue in water bridging) and Occupancy column (results for particular frame). Water molecules will be included in each frame.

In [33]: savePDBWaterBridges(waterBridges[0], coords_traj, PDBtraj_file[:-4]+'_frame0.pdb')

In [34]: savePDBWaterBridgesTrajectory(waterBridges, coords_traj,
   ....:                               filename=PDBtraj_file[:-4]+'_multi.pdb',
   ....:                               trajectory=trajectory)
   ....: 

Results saved in PDB file can be displayed as follows:

../../_images/Fig2.png