Cross-Docking and Scoring

[fisand08]

Dear Developers

As previously discussed I suggest some additions to your tool to make EasyDock even more easy.

It would be great to have some way of estimating pose prediction performance (i.e. by RMSD) in a re- or cross-docking scenario. A user could supply a set of X-ray structures with cocrystallized ligands and the best docking engine to reproduce them could be selected.

I typically spend a lot of time evaluating a particular algorithm, usually more than setting up the actual docking experiments to screen a library.

Best,
André

[Feriolet]

Hi André,

Thats interesting. Indeed it might be convenient to have some sort of way to perform cross/re docking easier. I did crossdocking for one of my projects in the past, and it took me a while to do it with easydock. Unfortunately, it's been a long time since I touch easydock github, and the main developers may have better idea.

If the logic of easydock has not been altered too much since I helped develop it, it may be difficult to accomodate this, as the logic was made to dock a set of ligands into a protein in one db file. Moreover, the repo mainly used the docking score instead of rmsd, so we may need to modify and add more parameter to the current ones. If you need to compare the rmsd between ligands across different proteins, imo I think it's better to use easydock as a python package, rather than the terminal one, since if i remember correctly, they have the mol_block or pdbqt string as an output.

Maybe Pavel have some better idea if we want to entertain this request?

[DrrDom]

Thank you for the suggestion,
indeed the current version does not automate these steps. The main feature of easydock - easy docking of a large set of compounds from command line and python using optionally a distributed environment. We extensively use easydock from python as a back-end in our automated design software crem-dock (https://github.com/ci-lab-cz/crem-dock) and crem-opt (not yet available) and therefore we try to keep this compatibility. However, what you requested can be automated.

I would like to more precisely specify use cases.

1. A user submit a set of ligands, a set of corresponding native proteins and reference ligand structures. As output it is expected to get RMSD for ligands? This would be a re-docking settings.
2. For cross-docking the input is the same, but the output will be RMSD for non-native proteins.

Both these cases can be run with the current version with a relatively simple shell script. If you have multiple directories where there are smi-file with a ligand, sdf/mol-file with the reference 3D structure of a ligand, a prepared protein, config.yml and grid box settings, you can run easydock using smi-file as input on a protein with a given grid box. If you specify --sdf argument it will return sdf-file with the top pose. All output sdf files from all directories can be combined in a single file. All reference ligand structures can also be combined in a single sdf file. As far as names of corresponding ligands in both files will be the same, you may use this script to calculate RMSD - https://github.com/DrrDom/rdkit-scripts/blob/master/rmsd_rdkit.py. The issue may be different protonation states in reference ligands and docked poses. In the ideal case, you may have reference ligands with correctly protonated centers and in this case docking can be run without protonation. Thus. ligands will perfectly match.
For cross-docking the script should be somewhat more complex, but the idea is the same.

The described solution may be a workaround. We used something very similar. As you may see at the presentation we did re-docking for 1400 complexes and it did not take much time. I think it should not be too difficult to implement this feature as a python script which will run on top of easydock. This seems the most reasonable approach, becasue we do not plan to change the database structure soon to adopt docking to multiple proteins (however, this is in a roadmap).

If somebody will implement this pipeline as a python script, I'll be happy to accept this contribution. I may look on this next week, but I cannot promise that we will do this quickly.

[GuzelMinibaeva]

Re-docking helper

directory layout

$BASE_DIR/
├─ target_001/
│ ├─ protein.pdbqt
│ ├─ config.yml
│ ├─ ligand.sdf
│ ├─ ligand.smi
│ └─ gridbox.txt
├─ target_002/
│ ├─ protein.pdbqt
│ ├─ config.yml
│ ├─ ligand.sdf
│ ├─ ligand.smi
│ └─ gridbox.txt

create_config.sh

#!/bin/bash

BASE_DIR="/ABSOLUTE/PATH/TO/BASE_DIR"

find "$BASE_DIR" -type d | while read -r dir; do
  protein_file=$(find "$dir" -maxdepth 1 -type f -name "*.pdbqt" ! -name "._*")  #There must be exactly one .pdbqt file for this to work
  gridbox_file="$dir/gridbox.txt"

  if [[ -f "${protein_file:-}" && -f "$gridbox_file" ]]; then
    config_file="$dir/config.yml"
    cat > "$config_file" <<EOF
protein: $protein_file
protein_setup: $gridbox_file
exhaustiveness: 16
seed: 120
n_poses: 1
ncpu: 4
EOF
    echo "Created: $config_file"
  else
    echo "Skipped: $dir"
  fi
done

run-docking.sh

#!/bin/bash

conda activate name_env

BASE_DIR="/ABSOLUTE/PATH/TO/BASE_DIR"

MAX_JOBS=4   #depends on available CPUs
JOBS=0

find "$BASE_DIR" -type d -links 2 | while read -r dir; do
    ligand_file=$(find "$dir" -maxdepth 1 -type f -name "*_ligand.smi" ! -name "._*" | head -n 1) #There must be exactly one .sdf or .smi file
    yml_file="$dir/vina_config.yml"
    db_file="${ligand_file%.smi}_dock.db"

    if [[ -f "$ligand_file" && -f "$yml_file" ]]; then
	
	run_dock -i "$ligand_file" -o "$db_file" --program vina --config "$yml_file" -c 1 --sdf &
	    
	JOBS=$((JOBS+1))

	if [ "$JOBS" -ge "$MAX_JOBS" ]; then
	    wait -n
	    JOBS=$((JOBS-1))
	fi
    fi
done

wait

combine SDFs to one file

cd $BASE_DIR

find . -maxdepth 2 -type f -name "*_dock.sdf" ! -name "._*" -size +0c -exec awk '1; END{if($0!="$$$$") print "$$$$"}' {} \; > ligands_from_complexes_docked.sdf