Restored bonds and clashed waters present after ligand preparation

[alissonws]

Hi!

I've been trying to prepare a .sdf guanosine molecule for a hydrated docking. However, after running mk_prepare_ligand.py, I can see a couple of water-water and water-ligand bonds created.

I tested this ligand in a docking run and it got worse. A lot of water molecules were either fused to the ligand, clashed with each other, or clashed inside residue chains.

Output ligand seems pretty much like this:

REMARK Flexibility Score: inf
ROOT
ATOM      1  C   UNL     1       3.386   0.023   0.423  1.00  0.00     0.178 C 
ATOM      2  O   UNL     1       2.161   0.738   0.655  1.00  0.00    -0.347 OA
ATOM      3  C   UNL     1       3.093  -1.054  -0.640  1.00  0.00     0.179 C 
ATOM      4  C   UNL     1       1.581  -0.897  -0.932  1.00  0.00     0.195 C 
ATOM      5  C   UNL     1       1.079  -0.164   0.338  1.00  0.00     0.253 C 
ATOM      6 WAT  UNL     1       2.054   2.990  -1.324  1.00  0.00     0.000 W 
ATOM      7 WAT  UNL     1       2.008   1.316   3.595  1.00  0.00     0.000 W 
ENDROOT
BRANCH   1   8
ATOM      8  C   UNL     1       4.460   0.983  -0.092  1.00  0.00     0.191 C 
BRANCH   8   9
ATOM      9  O   UNL     1       4.797   1.917   0.937  1.00  0.00    -0.394 OA
ATOM     10  H   UNL     1       5.473   2.558   0.680  1.00  0.00     0.210 HD
ATOM     11 WAT  UNL     1       6.010   0.293   3.149  1.00  0.00     0.000 W 
ATOM     12 WAT  UNL     1       2.307   3.513   1.439  1.00  0.00     0.000 W 
ATOM     13 WAT  UNL     1       6.872   3.885   0.148  1.00  0.00     0.000 W 
ENDBRANCH   8   9
ENDBRANCH   1   8
...

Any tips on why those bonds are being restored like this?

[rwxayheee]

Hi @alissonws, These are not actual bonds but water sites that will translate with the associated ligand atoms. Some visualization program might attempt to perceive them as W (Tungsten) atoms and try to connect the closely positioned W atoms. If you open it in PyMol, you won't see the connection between the water sites.

Once you get to the post-processing step after you use dry.py on your docking output, you will see which water sites may be replaced upon binding :)

[alissonws]

Oooh, that makes sense. So I shouldn't worry about these in my preparation files, right? One more thing, I'm also trying to run this protocol on flexible receptors. I've been following the steps and generating the maps and grid based on the rigid receptor, however, this leads to the misplacement of the ligand and the waters. Even after dry.py, there are ligand interactions with water molecules almost superposing flexible side chains, and the ligand is pushed to "unnatural" conformations. By "rigid" receptor, I mean the _rigid PDBQT file generated for flexible docking. The _flex part is used only when calling the vina binary. Any idea on what I'm missing?

Thanks!

[rwxayheee]

Hi @alissonws

there are ligand interactions with water molecules almost superposing flexible side chains

I think the water map 1uw6_receptor.W.map was computed without flexible sidechains (as if they were absent).
I don't know the best practice to resolve this.. Maybe you can regenerate the water map after docking, if you can merge the rigid receptor and the docked conformation of the flexible sidechains and run mapwater.py on the "combined" receptor. That way the water map should be aware of the flexible sidechains and not suggesting placing water right on these occupied spots.

[alissonws]

Hi, @rwxayheee , I agree. After a bit more reading, it seems we should manually remove WEAK waters if we find them misplaced. Even though only displaced waters are automatically removed on post-processing, clashing waters are scored down in the output file. This kind of clarifies the issue of generating a proper 3D model.

However, this post-processing will still be made over biased Vina results since the docking itself appears to miscalculate water-flexible-chain interactions. If the force field does not properly weigh down clashing waters, Vina results will be filled up with local minima poses that aren't even possible in real life. The docs itself says an experimentally validated pose is usually around 8 kcal/mol, so I believe the only solution is to increase --num_modes to a number we can properly screen around 8 kcal/mol poses while considering scores lower than that as artifacts of the scoring function. I assume we can reliably use the score value to compare poses after this manual curation. Does it make sense?

[rwxayheee]

Hi @alissonws,
I don't have enough experience with hydrated docking. I was suggesting to run dry.py with full-receptor pdbqt and maps, hoping that will do the removing ligand/receptor overlapping waters trick so the clashed waters can be removed by dry.py. I never did this before, maybe i should try it myself first..

I'm not very sure how the water-flexible-chain interaction is calculated in the docking calculations. Does W contribute to inter or intra energy? Can you help me understand why do you think Vina miscalculates it (did you get a weird inter or intra energy?