How to use AutoDock Vina as a component?

[adeptustao0bit]

I'm trying to use AutoDock Vina as a scoring component in REINVENT4, but it's not working as expected. Here's what I've done so far:

1. I placed the comp_vina.py file in the appropriate folder.
2. When I ran the command, I encountered a KeyError: 'vina' , which suggested that the component wasn't registered properly.
3. To troubleshoot, I used a script to help integrate the model. While the script runs without errors, all molecule affinity values generated by comp_vina are 0.00 .
4. I created an additional log file for comp_vina and noticed that the component isn't being called at all.
   I'm new to GitHub and would greatly appreciate any guidance on where the issue might be and how to resolve it. Thank you!

Here is my run_reinvent_with_plugin.py

import sys
import os
import subprocess

def setup_environment():
    """Set up the plugin environment and validate components"""
    # Add the plugin directory to the Python path
    plugin_base = r"D:\software\reinvent4.0\REINVENT4-main\reinvent_plugins\components"
    if plugin_base not in sys.path:
        sys.path.insert(0, plugin_base)

    # Validate the plugins
    try:
        import reinvent_plugins.components
        print("Components loaded successfully")

        from reinvent.scoring.importer import get_registry
        registry = get_registry()
        print(f"Found {len(registry)} components")

        # Check if the Vina component is registered
        vina_lookup = 'vina'.lower().replace("-", "").replace("_", "")
        if vina_lookup in registry:
            print(f"Vina registered: {vina_lookup}")
        else:
            print(f"Vina unregistered, searching for {vina_lookup}")
            print("Reachable components:", list(registry.keys()))

    except Exception as e:
        print(f"Plugin validation failed: {e}")
        return False

    return True

def run_reinvent(config_file, log_file=None):
    """Run the REINVENT4 framework with the scoring components"""

    if not setup_environment():
        print("Environment setup failed")
        return False

    # Build the command to run REINVENT4
    cmd = ["reinvent", config_file]
    if log_file:
        cmd.extend(["-l", log_file])

    print(f"Running command: {' '.join(cmd)}")

    try:
        # Execute the command
        process = subprocess.Popen(
            cmd,
            stdout=subprocess.PIPE,
            stderr=subprocess.STDOUT,  # Redirect stderr to stdout
            text=True,
            bufsize=1,  # Line buffering
            universal_newlines=True
        )

        # Capture and print output in real-time
        while True:
            output = process.stdout.readline()
            if output == '' and process.poll() is not None:
                break
            if output:
                print(output.strip())

        # Get the return code
        return_code = process.poll()

        if return_code == 0:
            print("REINVENT4 ran successfully")
            return True
        else:
            print(f"REINVENT4 failed with return code: {return_code}")
            return False

    except subprocess.CalledProcessError as e:
        print(f"REINVENT4 failed: {e}")
        print("STDERR:", e.stderr)
        return False


if __name__ == "__main__":
    # Replace these paths with your actual file paths
    config_file = r"D:\software\reinvent4.0\REINVENT4-main\configs\toml\staged_learning.toml"
    log_file = r"D:\software\reinvent4.0\REINVENT4-main\configs\toml\SL_generation\sl_1.log"

    run_reinvent(config_file, log_file)


#Here is my comp_vina.py

"""Run an external scoring subprocess for AutoDock Vina"""

from __future__ import annotations

__all__ = ["Vina"]

import os
import logging
import subprocess
import tempfile
from typing import List, Optional
from pathlib import Path

import numpy as np
from pydantic.dataclasses import dataclass
from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit.Chem.SaltRemover import SaltRemover
from rdkit.Chem.MolStandardize import rdMolStandardize
from rdkit.Chem.rdmolops import SanitizeFlags

from .component_results import ComponentResults
from .run_program import run_command
from .add_tag import add_tag
from reinvent_plugins.normalize import normalize_smiles

logger = logging.getLogger(__name__)


@add_tag("__parameters")
@dataclass
class Parameters:
    """Data class for storing parameters of the Vina scoring component"""
    acceptor_path: List[str]  # Path to the receptor structure file
    ligand_dir: List[str]  # Directory for temporary ligand files
    vina_path: List[str]  # Path to the AutoDock Vina executable
    log_dir: List[str]  # Directory for logging files
    output_dir: List[str]  # Directory for docking output files
    center_x: List[float]  # X-coordinate of the docking box center
    center_y: List[float]  # Y-coordinate of the docking box center
    center_z: List[float]  # Z-coordinate of the docking box center
    size_x: List[float]  # Size of the docking box along the X-axis
    size_y: List[float]  # Size of the docking box along the Y-axis
    size_z: List[float]  # Size of the docking box along the Z-axis


@add_tag("__component")
class Vina:
    """Class for running AutoDock Vina for molecular docking scoring"""

    def __init__(self, params: Parameters):
        # Initialize Vina with parameters
        self.vina_path = params.vina_path[0]
        self.ligand_dir = Path(params.ligand_dir[0])
        self.receptor_path = params.acceptor_path[0]
        self.log_dir = Path(params.log_dir[0])
        self.output_dir = Path(params.output_dir[0])

        # Docking box configuration
        self.center_x = params.center_x[0]
        self.center_y = params.center_y[0]
        self.center_z = params.center_z[0]
        self.size_x = params.size_x[0]
        self.size_y = params.size_y[0]
        self.size_z = params.size_z[0]

        # Create necessary directories
        self.ligand_dir.mkdir(parents=True, exist_ok=True)
        self.log_dir.mkdir(parents=True, exist_ok=True)
        self.output_dir.mkdir(parents=True, exist_ok=True)

        # Set up logging
        log_file_path = self.log_dir / "vina_plugin.log"
        self.logger = logging.getLogger('vina_plugin')
        self.logger.setLevel(logging.INFO)

        file_handler = logging.FileHandler(log_file_path)
        file_handler.setLevel(logging.INFO)

        formatter = logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s')
        file_handler.setFormatter(formatter)

        self.logger.addHandler(file_handler)

        self.logger.info("Vina component initialized successfully, log file created")

    @normalize_smiles
    def __call__(self, smiles_list: List[str]) -> ComponentResults:
        """Main scoring function called by REINVENT4 framework"""
        scores = []

        for idx, smiles in enumerate(smiles_list):
            try:
                score = self._process_single_molecule(smiles, idx)
                self.logger.info(f"Docking completed for SMILES: {smiles}")
                scores.append(score)
            except Exception as e:
                self.logger.error(f"Failed to process molecule {idx} with SMILES {smiles}: {e}")
                scores.append(0.0)  # Default score for failed molecules

        # Return results in the format required by REINVENT4
        return ComponentResults([np.array(scores, dtype=float)])

    def _process_single_molecule(self, smiles: str, molecule_id: int) -> float:
        """Process a single molecule through the Vina pipeline"""
        # Generate a unique filename for the molecule
        base_name = f"REINVENT_{self._get_internal_step()}_{molecule_id}"

        # Use a temporary directory to store intermediate files
        with tempfile.TemporaryDirectory() as temp_dir:
            temp_path = Path(temp_dir)

            # Convert SMILES to PDB format
            pdb_path = temp_path / f"{base_name}.pdb"
            if not self._smiles_to_pdb(smiles, str(pdb_path)):
                return 0.0
            self.logger.info(f"Converted SMILES to PDB format for molecule {smiles}")

            # Convert PDB to PDBQT format
            pdbqt_input_path = temp_path / f"{base_name}.pdbqt"
            if not self._convert_pdb_to_pdbqt(str(pdb_path), str(pdbqt_input_path)):
                return 0.0
            self.logger.info(f"Converted PDB to PDBQT format for molecule {smiles}")

            # Run AutoDock Vina docking
            output_path = self.output_dir / f"{base_name}_out.pdbqt"
            log_path = self.log_dir / f"{base_name}.log"

            affinity = self._run_vina_docking(
                str(pdbqt_input_path),
                str(output_path),
                str(log_path)
            )

            return affinity if affinity is not None else 0.0

    def _smiles_to_pdb(self, smiles: str, output_path: str, num_conformers: int = 50) -> bool:
        """Convert a SMILES string to a 3D PDB structure"""
        try:
            # Standardize the SMILES string
            standardized_smiles = self._standardize_smiles(smiles)
            if not standardized_smiles:
                return False

            # Create an RDKit molecule from the SMILES string
            mol = Chem.MolFromSmiles(standardized_smiles)
            if mol is None:
                return False

            # Add hydrogen atoms to the molecule
            mol = Chem.AddHs(mol)
            conformer_ids = AllChem.EmbedMultipleConfs(mol, numConfs=num_conformers, randomSeed=42)

            if not conformer_ids:
                return False

            # Optimize conformers and select the best one
            best_conf_id = self._optimize_and_select_conformer(mol, conformer_ids)
            if best_conf_id == -1:
                return False

            # Save the best conformer to a PDB file
            writer = Chem.PDBWriter(output_path)
            writer.write(mol, confId=best_conf_id)
            writer.close()

            return True

        except Exception as e:
            self.logger.error(f"Failed to convert SMILES to PDB: {e}")
            return False

    def _optimize_and_select_conformer(self, mol, conformer_ids, optimize_iterations: int = 100) -> int:
        """Optimize conformers and select the best one based on energy"""
        min_energy = float('inf')
        best_conf_id = -1

        for conf_id in conformer_ids:
            try:
                # Get molecule properties for MMFF optimization
                mmff_props = AllChem.MMFFGetMoleculeProperties(mol, mmffVariant='MMFF94s')
                if mmff_props is None:
                    continue

                # Optimize the conformer
                result = AllChem.MMFFOptimizeMoleculeConfs(
                    mol,
                    maxIters=optimize_iterations,
                    mmffVariant='MMFF94s'
                )

                # Update the best conformer if a lower energy is found
                for i, (not_converged, energy) in enumerate(result):
                    if not_converged == 0 and energy < min_energy:
                        min_energy = energy
                        best_conf_id = i

            except Exception as e:
                logger.debug(f"Optimization failed for conformer {conf_id}: {e}")
                continue

        return best_conf_id

    def _convert_pdb_to_pdbqt(self, pdb_path: str, pdbqt_path: str) -> bool:
        """Convert a PDB file to PDBQT format using Open Babel"""
        try:
            # Use Open Babel to perform the conversion
            command = [
                "obabel", pdb_path, "-O", pdbqt_path,
                "--partialcharge", "gasteiger"
            ]

            subprocess.run(
                command,
                check=True,
                capture_output=True,
                text=True
            )

            # Check if the output file was created successfully
            return os.path.exists(pdbqt_path)

        except subprocess.CalledProcessError as e:
            self.logger.error(f"Failed to convert PDB to PDBQT: {e}")
            return False
        except FileNotFoundError:
            self.logger.error("Open Babel (obabel) not found in PATH")
            return False

    def _run_vina_docking(self, ligand_path: str, output_path: str, log_path: str) -> Optional[float]:
        """Run AutoDock Vina docking and parse the results"""
        try:
            # Prepare the command for AutoDock Vina
            command = [
                self.vina_path,
                "--receptor", self.receptor_path,
                "--ligand", ligand_path,
                "--center_x", str(self.center_x),
                "--center_y", str(self.center_y),
                "--center_z", str(self.center_z),
                "--size_x", str(self.size_x),
                "--size_y", str(self.size_y),
                "--size_z", str(self.size_z),
                "--out", output_path,
                "--log", log_path,
                "--seed", "42"
            ]

            # Execute the AutoDock Vina command
            run_command(command)

            # Parse the docking results
            return self._parse_vina_log(log_path)

        except Exception as e:
            logger.error(f"Failed to run AutoDock Vina: {e}")
            return None

    def _parse_vina_log(self, log_path: str) -> Optional[float]:
        """Parse the Vina log file to extract the binding affinity"""
        try:
            # Open and read the log file
            with open(log_path, 'r') as log_file:
                for line in log_file:
                    if line.strip().startswith('1 '):  # The first line contains the best binding mode
                        parts = line.split()
                        if len(parts) >= 2:
                            affinity = float(parts[1])
                            logger.info(f"Found binding energy: {affinity} kcal/mol")
                            return affinity
            self.logger.warning(f"No binding energy found in log file: {log_path}")
            return None
        except (IOError, ValueError, IndexError) as e:
            self.logger.error(f"Failed to parse Vina log file: {e}")
            return None

    def _standardize_smiles(self, smiles: str) -> Optional[str]:
        """Standardize a SMILES string using RDKit"""
        try:
            # Create an RDKit molecule from the SMILES string
            mol = Chem.MolFromSmiles(smiles)
            if mol is None:
                return None

            # Remove salts from the molecule
            salt_remover = SaltRemover()
            desalted_mol = salt_remover.StripMol(mol, dontRemoveEverything=True)

            # Validate the molecule
            if not self._validate_molecule(desalted_mol):
                return None

            # Sanitize the molecule
            Chem.SanitizeMol(
                desalted_mol,
                SanitizeFlags.SANITIZE_ALL ^ SanitizeFlags.SANITIZE_CLEANUPCHIRALITY
            )

            # Normalize the molecule
            normalizer = rdMolStandardize.Normalizer()
            standardized_mol = normalizer.normalize(desalted_mol)

            # Return the standardized SMILES string
            return Chem.MolToSmiles(
                standardized_mol,
                canonical=True,
                isomericSmiles=True,
                kekuleSmiles=True
            )

        except Exception as e:
            self.logger.error(f"Failed to standardize SMILES: {e}")
            return None

    def _validate_molecule(self, mol: Chem.Mol) -> bool:
        """Check if the molecule is valid"""
        return mol.GetNumAtoms() > 0 and any(atom.GetAtomicNum() > 0 for atom in mol.GetAtoms())

    def _get_internal_step(self) -> int:
        """Generate a unique identifier for file naming"""
        # Use a timestamp to ensure uniqueness
        import time
        return int(time.time() * 1000)

#Here is the soring part of staged_learning.toml
[stage.scoring]
type = "geometric_mean"  # aggregation function

[[stage.scoring.component]]
# Custom alerts if used in a custom_product filter out unwanted groups
[stage.scoring.component.custom_alerts]

[[stage.scoring.component.custom_alerts.endpoint]]
name = "Unwanted SMARTS"  # user chosen name for output
weight = 0.79  # weight to fine-tune the relevance of this component

# parameters for the component:
# a list of unwanted SMARTS(!) to be scored as zero
params.smarts = [
    "[*;r8]",
    "[*;r9]",
    "[*;r10]",
    "[*;r11]",
    "[*;r12]",
    "[*;r13]",
    "[*;r14]",
    "[*;r15]",
    "[*;r16]",
    "[*;r17]",
    "[#8][#8]",
    "[#6;+]",
    "[#16][#16]",
    "[#7;!n][S;!$(S(=O)=O)]",
    "[#7;!n][#7;!n]",
    "C#C",
    "C(=[O,S])[O,S]",
    "[#7;!n][C;!$(C(=[O,N])[N,O])][#16;!s]",
    "[#7;!n][C;!$(C(=[O,N])[N,O])][#7;!n]",
    "[#7;!n][C;!$(C(=[O,N])[N,O])][#8;!o]",
    "[#8;!o][C;!$(C(=[O,N])[N,O])][#16;!s]",
    "[#8;!o][C;!$(C(=[O,N])[N,O])][#8;!o]",
    "[#16;!s][C;!$(C(=[O,N])[N,O])][#16;!s]"
]

[[stage.scoring.component]]
[stage.scoring.component.MolecularWeight]

[[stage.scoring.component.MolecularWeight.endpoint]]
name = "Molecular weight"  # user chosen name for output
weight = 0.342  # weight to fine-tune the relevance of this component

# A transform ensures that the output from the scoring component ranges
# from 0 to 1 to serve as a proper score.  Here we use a double sigmoid
# to transform weights into the range 200-500 a.u.
transform.type = "double_sigmoid"
transform.high = 500.0
transform.low = 200.0
transform.coef_div = 500.0
transform.coef_si = 20.0
transform.coef_se = 20.0

[[stage.scoring.component]]
[stage.scoring.component.Vina]
[[stage.scoring.component.Vina.endpoint]]
name = "Affinity"
weight = 0.80

params.acceptor_path = "D:/study/python/autodock/acceptor/acceptor.pdbqt"
params.ligand_dir = "D:/software/reinvent4.0/REINVENT4-main/configs/toml/SL_generation/ligand"
params.vina_path = "D:/software/reinvent4.0/REINVENT4-main/configs/toml/SL_genertaion/vina.exe"
params.log_dir = "D:/software/reinvent4.0/REINVENT4-main/configs/toml/SL_generation/log"
params.output_dir = "D:/software/reinvent4.0/REINVENT4-main/configs/toml/SL_generation/output"
params.center_x = 36.6
params.center_y = -70
params.center_z = 6.9
params.size_x = 50
params.size_y = 50
params.size_z = 50

[halx]

Hi,

many thanks for the detailed report. First, make sure that you do not have a init.py file in the directory with your component. Second, have you had a look into maize which comes with various docking workflows?

Many thanks,
Hannes.

[adeptustao0bit]

really thanks for your fast reply, I will have a look into maize.