bio-applied-virtual-screening

name	bio-applied-virtual-screening
description	Dock ligand libraries with AutoDock Vina/meeko, filter hits by ADMET (Lipinski, LogS, hERG), rank by composite docking+QSAR score in pandas. Use when docking SMILES/SDF vs a target or prioritizing virtual screening hits.
tool_type	python
primary_tool	AutoDock Vina

Virtual Screening & ADMET Prediction

When to Use

• Screening a compound library (SMILES/SDF) against a protein target by molecular docking
• Preparing a receptor (PDB) and ligands (SDF) for AutoDock Vina
• Filtering docking hits with ADMET rules (Lipinski Ro5, LogS, CYP inhibition, hERG risk)
• Combining docking score + QSAR prediction + ADMET pass/fail into one ranked hit list
• Building a structure-based virtual screening pipeline for hit-to-lead prioritization

Version Compatibility

AutoDock Vina ≥1.2, meeko ≥0.5 (ligand/receptor PDBQT prep), RDKit ≥2023.09, pandas ≥2.0, Python ≥3.10. Vina 1.2+ takes PDBQT for both receptor and ligand (no separate prepare_receptor4.py/prepare_ligand4.py MGLTools scripts needed).

Prerequisites

• pip install vina meeko rdkit pandas (or conda install -c conda-forge vina)
• A cleaned receptor structure (waters/ligands stripped) and a defined binding-site box
• Familiarity with bio-chemoinformatics-molecular-io (SMILES/SDF handling) and bio-structural-biology-structure-io (PDB parsing)

Protein Structure Preparation

Goal: get a receptor PDBQT with hydrogens added and no waters/crystallographic ligands.

Approach: download the PDB, strip heteroatoms, add hydrogens with PDBFixer, then convert to PDBQT with meeko/OpenBabel. The classic Vina tutorial target is PDB 1IEP — the Abl kinase (ABL1) tyrosine kinase domain in complex with imatinib — not EGFR.

# Download and fix the Abl kinase-imatinib structure (PDB 1IEP)
wget -q https://files.rcsb.org/download/1IEP.pdb
pdbfixer 1IEP.pdb --output 1IEP_fixed.pdb --add-hydrogens --remove-heterogens --keep-water=False

# Convert receptor PDB to PDBQT (Vina 1.2+ / meeko)
mk_prepare_receptor.py -i 1IEP_fixed.pdb -o 1IEP_fixed --box_size 20 20 20 \
    --box_center 22.5 5.0 18.0

from Bio.PDB import PDBParser, PDBIO, Select


class HetatmWaterSelect(Select):
    """Keep only standard amino-acid residues; drop waters and heteroatoms."""

    def accept_residue(self, residue):
        return residue.id[0] == " "  # blank hetfield = standard residue


def strip_to_protein_only(pdb_in: str, pdb_out: str) -> None:
    """Write a protein-only PDB, removing waters/ligands/ions before docking prep."""
    parser = PDBParser(QUIET=True)
    structure = parser.get_structure("receptor", pdb_in)
    io = PDBIO()
    io.set_structure(structure)
    io.save(pdb_out, select=HetatmWaterSelect())

Molecular Docking with AutoDock Vina

Goal: dock a ligand (or library of ligands) into the prepared binding site and get a binding-affinity score (kcal/mol; more negative = better predicted binding).

Approach: convert each ligand SDF/SMILES to PDBQT with meeko, then run Vina against the receptor with a box centered on the known active site.

# Ligand prep: SDF -> PDBQT (adds Gasteiger charges, sets rotatable bonds)
mk_prepare_ligand.py -i imatinib.sdf -o imatinib.pdbqt

# Run Vina docking against the Abl kinase ATP-binding site
vina --receptor 1IEP_fixed.pdbqt --ligand imatinib.pdbqt \
     --center_x 22.5 --center_y 5.0 --center_z 18.0 \
     --size_x 20 --size_y 20 --size_z 20 \
     --exhaustiveness 8 \
     --out imatinib_docked.pdbqt --log docking.log

grep "REMARK VINA RESULT" imatinib_docked.pdbqt | head -1

from vina import Vina


def dock_ligand(receptor_pdbqt: str, ligand_pdbqt: str, center: tuple, box_size=(20, 20, 20)) -> float:
    """Dock one ligand with the Vina Python API; return the best (lowest) binding score in kcal/mol."""
    v = Vina(sf_name="vina")
    v.set_receptor(receptor_pdbqt)
    v.set_ligand_from_file(ligand_pdbqt)
    v.compute_vina_maps(center=list(center), box_size=list(box_size))
    v.dock(exhaustiveness=8, n_poses=10)
    energies = v.energies(n_poses=1)  # [[total, inter, intra, torsion, intra_best]]
    return float(energies[0][0])

ADMET Property Prediction

Goal: flag compounds with poor drug-likeness or predicted toxicity before/after docking.

Approach: compute cheap RDKit-based Lipinski Ro5 descriptors as a first-pass filter; use SwissADME (web) or a DeepChem Tox21 model for hERG/CYP/toxicity predictions on the survivors.

from rdkit import Chem
from rdkit.Chem import Descriptors, Lipinski


def passes_lipinski(smiles: str, max_violations: int = 1) -> bool:
    """Return True if a SMILES passes Lipinski's Rule of Five (allowing max_violations)."""
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return False
    violations = sum([
        Descriptors.MolWt(mol) > 500,
        Descriptors.MolLogP(mol) > 5,
        Lipinski.NumHDonors(mol) > 5,
        Lipinski.NumHAcceptors(mol) > 10,
    ])
    return violations <= max_violations

Hit Prioritization

Goal: merge docking score, QSAR-predicted pIC50, and ADMET pass/fail into one ranked shortlist.

Approach: build a pandas DataFrame per compound, filter on admet_pass, then rank by a composite score (flip the sign on vina_score since more negative = better binding, so higher composite = better hit).

import pandas as pd


def prioritize_hits(smiles_library, docking_scores, qsar_predictions, admet_flags, top_n=10) -> pd.DataFrame:
    """Combine Vina score + QSAR pIC50 + ADMET pass/fail into a ranked top-N hit table."""
    hits = pd.DataFrame({
        "SMILES": smiles_library,
        "vina_score": docking_scores,      # kcal/mol, more negative = stronger predicted binding
        "qsar_pIC50": qsar_predictions,     # higher = more potent
        "admet_pass": admet_flags,         # bool: passed ADMET filter
    })
    hits_filtered = hits[hits["admet_pass"]].copy()
    hits_filtered["composite"] = -hits_filtered["vina_score"] + hits_filtered["qsar_pIC50"]
    top_hits = hits_filtered.nlargest(top_n, "composite")
    return top_hits[["SMILES", "vina_score", "qsar_pIC50", "composite"]]

Pitfalls

• Docking score ≠ binding affinity: Vina scores correlate weakly with experimental Kd/IC50; always validate top hits experimentally before further investment
• Receptor flexibility: rigid docking misses induced-fit binding; use ensemble docking (multiple receptor conformations) for flexible targets
• PAINS compounds: filter pan-assay interference structures (rdkit.Chem.FilterCatalog) before scoring — they show up as false-positive hits across unrelated assays
• Box placement: an incorrectly centered/sized docking box silently produces poses outside the real binding site with plausible-looking scores
• Protonation state: skipping hydrogen addition/pH-appropriate tautomers on both receptor and ligand biases scores and pose geometry

See Also

• bio-chemoinformatics-virtual-screening — ligand-based screening (similarity, pharmacophore, QSAR) as an alternative/complement to docking
• bio-chemoinformatics-admet-prediction — deeper ADMET modeling (LogS, CYP, hERG) beyond Lipinski filters
• bio-structural-biology-structure-io — general PDB/mmCIF parsing and structure cleanup
• bio-chemoinformatics-molecular-io — SMILES/SDF reading, writing, and format conversion