<p>I have the results of a virtual screening experiment using docking simulations with Autodock Vina. The result is a matrix of 7 (proteins) by 28000 (ligands) with the calculated binding energies for each protein/ligand pair:</p>
<div class="s-table-container">
<table class="s-table">
<thead>
<tr>
<th>Ligands - Proteins</th>
<th>Protein1</th>
<th>Protein2</th>
<th>Protein3</th>
<th>Protein4</th>
<th>Protein5</th>
<th>Protein6</th>
<th>Protein7</th>
</tr>
</thead>
<tbody>
<tr>
<td>Ligand1</td>
<td>-5.4</td>
<td>-5.9</td>
<td>-5.0</td>
<td>-5.2</td>
<td>-5.8</td>
<td>-5.7</td>
<td>-4.5</td>
</tr>
<tr>
<td>Ligand2</td>
<td>-9.5</td>
<td>-8.6</td>
<td>-8.9</td>
<td>-9.0</td>
<td>-9.2</td>
<td>-7.1</td>
<td>-6.4</td>
</tr>
<tr>
<td>Ligand3</td>
<td>-9.5</td>
<td>-8.3</td>
<td>-8.2</td>
<td>-9.9</td>
<td>nan</td>
<td>-6.8</td>
<td>nan</td>
</tr>
<tr>
<td>Ligand4</td>
<td>-9.4</td>
<td>-7.8</td>
<td>-8.0</td>
<td>-6.9</td>
<td>-9.5</td>
<td>-1.3</td>
<td>-7.8</td>
</tr>
<tr>
<td>Ligand5</td>
<td>-9.3</td>
<td>-6.9</td>
<td>-8.5</td>
<td>-7.9</td>
<td>-9.2</td>
<td>-6.8</td>
<td>-6.9</td>
</tr>
<tr>
<td>---</td>
<td></td>
<td></td>
<td></td>
<td></td>
<td></td>
<td></td>
<td></td>
</tr>
<tr>
<td>Ligand28000</td>
<td>-6.3</td>
<td>-7.9</td>
<td>-5.4</td>
<td>-2.7</td>
<td>-10.1</td>
<td>-9.4</td>
<td>-3.9</td>
</tr>
</tbody>
</table>
</div>
<p>nan results are either failed dockings, or positive values.</p>
<p>The scores are not only affected by how well a ligand docks with each protein, but also <strong>bigger (higher molecular weight) ligands tend to show better (more negative) scores</strong>. However, this is in part an artifact and it does not imply that bigger molecules are better and smaller ones are worse.</p>
<p><a href="https://i.sstatic.net/JgNKo.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/JgNKo.png" alt="energy/size graph" /></a></p>
<p>Typically, for a virtual screening with just one protein, you would select the top ranking ligands and that would be it (e.g. keep the top 1% and test them). However, this results in some ligands selected for most proteins because they invariably score better due to their size.</p>
<p>In this case, I could also take into account the variability of each ligand because I have observations for different unrelated proteins (i.e. see if a particular ligand systematically binds with better scores across proteins).</p>
<p><strong>Is there any way to take into consideration both factors when selecting the best ligands?</strong></p>
<p>I thought of calculating Z-scores 1) across rows and 2) across columns, and select the ligands for each protein that have less than -1.65 standard deviations (p=0.05) and that also fill the same criteria across proteins. However, there are only 7 observations for each ligand to calculate the Z-scores, and it does not feel correct to give it the same weight to the Z-scores calculated from the 28000 observations for each protein. Plus, in the end, protein docking effect is more important than the ligand effect.</p>
<p>Any suggestions would be greatly appreciated!</p>
<p><a href="https://jcheminf.biomedcentral.com/articles/10.1186/s13321-019-0330-2" rel="nofollow noreferrer">Ligand efficiency</a> looks like a promising suggestion. However, upon trying (calculating ∆G/MW or ∆G/n.heavy_atoms) resulted in very small ligands with mediocre or bad energies going to the top (e.g. urea with -3.1 became a Z-score of -8.5 SDs).</p>

CC BY-SA 4.0
2022-09-02T10:37:58.407Z
Urea?! —what dataset are you docking and what's the application? I had assumed it was Enamine REAL or Wuxi biologics dataset for purchasing for drug discovery.
Matteo Ferla
CC BY-SA 4.0
2022-09-02T10:44:55.190Z
If it is to find the native ligand of an enzyme, the task is completely different, but the monetary aspect is worse —say the native ligand is something like phosphoribosyl pyrophosphate and you have 10 candidates that are close, that will cost as much as the GDP of a small state. Therefore ancillary checks are needed like operon guilty-by-association etc.
Matteo Ferla
CC BY-SA 4.0
2022-09-02T10:46:45.327Z
I seem to have glitched the system by adding an extra comment and my previous comment has shrunk. If it's for drug discovery I'll re type it... :shrung:
Matteo Ferla
CC BY-SA 4.0
2022-09-02T11:04:21.573Z
Sorry Matteo I don't know where my comment went either! I used the ZINC database of natural compounds [here](https://zinc.docking.org/substances/subsets/natural-products/), with the aim of drug discovery. It is expected to find things like urea I guess, and of course any result like this would be manually curated
albertr
CC BY-SA 4.0
2022-09-02T11:40:23.400Z
My comment was about cost and success rate: smaller compounds are cheaper and more likely to be true so LE + a lower cutoff of 150–200 Da depending on the assay or co-crystal soak. So generally smaller compound hits are elaborated upon until one has a lead and jumping to the final product is a waste of cash. ZINC natural products is not make-on-demand & may require a chemist to source the plant (read: import paperwork) and extract the product (read: low yield). Plus they may not be not "drug-like" (cf. Lipinski's Ro5 vs. nature's proclivity towards PO4, CoA or glycosyl groups as a handle).
Matteo Ferla
CC BY-SA 4.0
2022-09-02T11:53:27.443Z
What I did was to keep any ZINC "natural-products" entry also described as "for-sale", which (in theory...) should be available for buying from distributors, and also discard any that violated 2 or more Lipinski's rules. But yeah, still it left with around 28000 compounds. I will try to see if I can use LE without boosting the "very" small compounds.
albertr
CC BY-SA 4.0
2022-09-02T13:34:36.833Z
Ah. 28k is more than I thought. I find natural compounds awesome to stare at to see their biosynthesis, but I'd worry the routes for expansions may be more isolated than the modular building block approach —say you have a cyclised sesquiterpene (farnesol (3x isoprene) is 222 Da) and you would like to change a substituent it would need to exist in nature or on sale or easily modifiable (Codexis & co. sell plates of exciting enzymes but AFIK that is for production optimisation not exploration).
Matteo Ferla
<r><docking><drugs><cheminformatics>
How to select the best ligands in a virtual screening matrix
# How to select the best ligands in a virtual screening matrix

(Asked by: username_0 on 2022-09-01)

I have the results of a virtual screening experiment using docking simulations with Autodock Vina. The result is a matrix of 7 (proteins) by 28000 (ligands) with the calculated binding energies for each protein/ligand pair:

| Ligands - Proteins | Protein1 | Protein2 | Protein3 | Protein4 | Protein5 | Protein6 | Protein7 |
| --- | --- | --- | --- | --- | --- | --- | --- |
| Ligand1 | -5.4 | -5.9 | -5.0 | -5.2 | -5.8 | -5.7 | -4.5 |
| Ligand2 | -9.5 | -8.6 | -8.9 | -9.0 | -9.2 | -7.1 | -6.4 |
| Ligand3 | -9.5 | -8.3 | -8.2 | -9.9 | nan | -6.8 | nan |
| Ligand4 | -9.4 | -7.8 | -8.0 | -6.9 | -9.5 | -1.3 | -7.8 |
| Ligand5 | -9.3 | -6.9 | -8.5 | -7.9 | -9.2 | -6.8 | -6.9 |
| --- |  |  |  |  |  |  |  |
| Ligand28000 | -6.3 | -7.9 | -5.4 | -2.7 | -10.1 | -9.4 | -3.9 |

nan results are either failed dockings, or positive values.

The scores are not only affected by how well a ligand docks with each protein, but also **bigger (higher molecular weight) ligands tend to show better (more negative) scores**. However, this is in part an artifact and it does not imply that bigger molecules are better and smaller ones are worse.

[![energy/size graph](https://i.sstatic.net/JgNKo.png)](https://i.sstatic.net/JgNKo.png)

Typically, for a virtual screening with just one protein, you would select the top ranking ligands and that would be it (e.g. keep the top 1% and test them). However, this results in some ligands selected for most proteins because they invariably score better due to their size.

In this case, I could also take into account the variability of each ligand because I have observations for different unrelated proteins (i.e. see if a particular ligand systematically binds with better scores across proteins).

**Is there any way to take into consideration both factors when selecting the best ligands?**

I thought of calculating Z-scores 1) across rows and 2) across columns, and select the ligands for each protein that have less than -1.65 standard deviations (p=0.05) and that also fill the same criteria across proteins. However, there are only 7 observations for each ligand to calculate the Z-scores, and it does not feel correct to give it the same weight to the Z-scores calculated from the 28000 observations for each protein. Plus, in the end, protein docking effect is more important than the ligand effect.

Any suggestions would be greatly appreciated!

[Ligand efficiency](https://jcheminf.biomedcentral.com/articles/10.1186/s13321-019-0330-2) looks like a promising suggestion. However, upon trying (calculating ∆G/MW or ∆G/n.heavy\_atoms) resulted in very small ligands with mediocre or bad energies going to the top (e.g. urea with -3.1 became a Z-score of -8.5 SDs).

**Tags:** r, docking, drugs, cheminformatics
> **Comments:**
> **username_1**: Urea?! —what dataset are you docking and what's the application? I had assumed it was Enamine REAL or Wuxi biologics dataset for purchasing for drug discovery.  
> **username_1**: If it is to find the native ligand of an enzyme, the task is completely different, but the monetary aspect is worse —say the native ligand is something like phosphoribosyl pyrophosphate and you have 10 candidates that are close, that will cost as much as the GDP of a small state. Therefore ancillary checks are needed like operon guilty-by-association etc.  
> **username_1**: I seem to have glitched the system by adding an extra comment and my previous comment has shrunk. If it's for drug discovery I'll re type it... :shrung:  
> **username_0**: Sorry Matteo I don't know where my comment went either! I used the ZINC database of natural compounds [here](https://zinc.docking.org/substances/subsets/natural-products/), with the aim of drug discovery. It is expected to find things like urea I guess, and of course any result like this would be manually curated  
> **username_1**: My comment was about cost and success rate: smaller compounds are cheaper and more likely to be true so LE + a lower cutoff of 150–200 Da depending on the assay or co-crystal soak. So generally smaller compound hits are elaborated upon until one has a lead and jumping to the final product is a waste of cash. ZINC natural products is not make-on-demand & may require a chemist to source the plant (read: import paperwork) and extract the product (read: low yield). Plus they may not be not "drug-like" (cf. Lipinski's Ro5 vs. nature's proclivity towards PO4, CoA or glycosyl groups as a handle).  
> **username_0**: What I did was to keep any ZINC "natural-products" entry also described as "for-sale", which (in theory...) should be available for buying from distributors, and also discard any that violated 2 or more Lipinski's rules. But yeah, still it left with around 28000 compounds. I will try to see if I can use LE without boosting the "very" small compounds.  
> **username_1**: Ah. 28k is more than I thought. I find natural compounds awesome to stare at to see their biosynthesis, but I'd worry the routes for expansions may be more isolated than the modular building block approach —say you have a cyclised sesquiterpene (farnesol (3x isoprene) is 222 Da) and you would like to change a substituent it would need to exist in nature or on sale or easily modifiable (Codexis & co. sell plates of exciting enzymes but AFIK that is for production optimisation not exploration).  

---

## 0 answers

