Lesson 11 of 12 · 11 min

An Introduction to Virtual Screening

So far you have docked a single ligand into one binding site. Real drug-discovery projects flip that around and push a whole library of candidate compounds through the same receptor, then rank them by predicted affinity. That process is virtual screening, and its job is triage: cheaply narrowing thousands of molecules down to a shortlist worth studying further.

Prepare the library

bash

mkdir -p pdbqt
# ligands.smi holds one "SMILES<TAB>name" per line
while read -r smiles name; do
  printf '%s\n' "$smiles" | obabel -ismi -O pdbqt/"$name".pdbqt --gen3d -p 7.4
done < ligands.smi

Tip: OpenBabel assigns a single protonation state and one tautomer per molecule. The -p 7.4 flag adds hydrogens for physiological pH, but it is only an approximation. For a serious campaign, enumerate tautomers and protomers with a dedicated tool before docking, since binding can hinge on the exact charge state.

Batch-dock with Vina

The receptor and search box do not change across the library, so define them once in a config file and reuse it for every ligand. Loop Vina over each PDBQT, writing one pose file and one log per compound. Screening is throughput-bound, so keep exhaustiveness modest and dock several compounds in parallel when you have the cores.

bash

cat > conf.txt <<'EOF'
receptor = receptor.pdbqt
center_x = 11.0
center_y = 92.5
center_z = 24.0
size_x = 22
size_y = 22
size_z = 22
exhaustiveness = 8
EOF

mkdir -p out logs
for lig in pdbqt/*.pdbqt; do
  name=$(basename "$lig" .pdbqt)
  vina --config conf.txt --ligand "$lig" \
       --out out/"${name}_out.pdbqt" > logs/"$name".log 2>&1
done

bash

printf "compound\taffinity\n" > ranked.tsv
for pose in out/*_out.pdbqt; do
  name=$(basename "$pose" _out.pdbqt)
  best=$(grep -m1 "REMARK VINA RESULT" "$pose" | awk '{print $4}')
  printf "%s\t%s\n" "$name" "$best"
done | sort -k2,2n >> ranked.tsv

column -t ranked.tsv

compound      affinity
compound_042  -10.3
compound_311  -9.8
compound_007  -9.5
compound_128  -9.1
compound_256  -8.7
compound_019  -8.2
compound_204  -7.6
compound_088  -6.9

The best pose per compound is its mode-1 affinity in kcal/mol, and more negative means a stronger predicted binder. Sorting ascending puts your enriched set at the top; a common first cut keeps everything below a threshold such as -8 kcal/mol, or simply the top few percent. Treat that shortlist as hypotheses, not confirmed hits.

Warning: docking scores are noisy. Vina's scoring function omits explicit water, treats the receptor as rigid, and only crudely approximates entropy (mostly a rotatable-bond penalty), so it produces many false positives and the ranking is unreliable near the cutoff. Use screening only to prioritize, then confirm survivors by visually inspecting the poses, rescoring or running short MD, and ultimately wet-lab assays.

Try it yourself: Grab a small set of 10 to 20 ligands as SMILES (for example a ZINC or Enamine subset), convert them to PDBQT with the OpenBabel loop above, and batch-dock them into a receptor you have already prepared. Build ranked.tsv, then report your top three compounds by affinity and open the best pose in PyMOL to check that it actually sits inside the pocket.

Visualizing Docking Results in PyMOL Wrap-Up and Where to Go Next