Protein Structure Prediction for Drug Discovery

AlphaFold2 (2021) fundamentally changed structure-based drug discovery by making accurate protein structure prediction accessible for any protein sequence. AlphaFold3 (2024) extended this to protein-ligand, protein-nucleic acid, and protein-protein complexes. ESMFold provides single-sequence prediction (no MSA required) with lower accuracy but ~60× faster inference — crucial for rapid screening of protein variant panels or novel targets.

The critical caveat for drug discovery: structure prediction accuracy varies dramatically by region. Using AlphaFold structures without quality-gating by pLDDT leads to docking artefacts and false positive hits. This guide covers rigorous quality filtering and docking preparation for predicted structures.

pLDDT and PAE: Quality Metrics

AlphaFold provides two per-residue/pair confidence metrics:

• pLDDT (predicted Local Distance Difference Test): per-residue confidence 0–100. >90 = very high confidence (comparable to crystal structure). 70–90 = good. 50–70 = low confidence (treat with caution). <50 = disordered region.
• PAE (Predicted Aligned Error): for protein complexes (AlphaFold3/Multimer), the expected positional error in Å for residue i when residue j is used as the alignment reference. Low inter-domain PAE indicates reliable domain-domain orientation.

import numpy as np
import json
from pathlib import Path

def parse_alphafold_json(af_json_path: str) -> dict:
    """
    Parse AlphaFold2 result JSON to extract pLDDT scores and PAE matrix.

    Args:
        af_json_path: Path to AlphaFold result_model_1.json (or .pkl unpickled)

    Returns:
        Dict with per-residue pLDDT and PAE matrix (if present).
    """
    with open(af_json_path) as f:
        data = json.load(f)

    plddt = np.array(data.get("plddt", []))
    pae = np.array(data.get("predicted_aligned_error", []))

    return {
        "plddt": plddt,
        "pae": pae,
        "mean_plddt": float(plddt.mean()) if len(plddt) > 0 else None,
        "n_residues": len(plddt),
    }


def assess_structure_quality(plddt: np.ndarray) -> dict:
    """
    Classify structure quality based on per-residue pLDDT values.

    Returns:
        Dict with residue counts per confidence tier and druggability verdict.
    """
    very_high = int((plddt >= 90).sum())
    high = int(((plddt >= 70) & (plddt < 90)).sum())
    low = int(((plddt >= 50) & (plddt < 70)).sum())
    disordered = int((plddt < 50).sum())
    total = len(plddt)

    pct_well_defined = (very_high + high) / total * 100 if total > 0 else 0

    verdict = (
        "suitable_for_docking" if pct_well_defined >= 70
        else "caution_significant_disordered_regions" if pct_well_defined >= 40
        else "not_recommended_for_structure_based_design"
    )

    return {
        "n_very_high_plddt_ge90": very_high,
        "n_high_plddt_70_90": high,
        "n_low_plddt_50_70": low,
        "n_disordered_plddt_lt50": disordered,
        "pct_well_defined": round(pct_well_defined, 1),
        "verdict": verdict,
    }


def get_reliable_binding_site_residues(
    plddt: np.ndarray,
    binding_site_residue_indices: list[int],
    min_plddt: float = 70.0,
) -> dict:
    """
    Check whether the predicted binding site residues meet quality threshold.

    Args:
        plddt: Per-residue pLDDT array (0-indexed)
        binding_site_residue_indices: 0-indexed residue positions forming the pocket
        min_plddt: Minimum acceptable pLDDT for binding site residues

    Returns:
        Quality assessment of binding site residues.
    """
    site_plddt = plddt[binding_site_residue_indices]
    n_reliable = int((site_plddt >= min_plddt).sum())
    pct_reliable = n_reliable / len(site_plddt) * 100 if len(site_plddt) > 0 else 0

    return {
        "site_residue_count": len(binding_site_residue_indices),
        "n_reliable_site_residues": n_reliable,
        "pct_reliable": round(pct_reliable, 1),
        "mean_site_plddt": round(float(site_plddt.mean()), 1),
        "min_site_plddt": round(float(site_plddt.min()), 1),
        "site_reliable": pct_reliable >= 80,
    }

Preparing Predicted Structures for Docking

from Bio.PDB import PDBParser, PDBIO, Select
import numpy as np

class HighConfidenceSelect(Select):
    """BioPython Select subclass to filter residues below pLDDT threshold."""

    def __init__(self, plddt_array: np.ndarray, min_plddt: float = 50.0):
        self.plddt = plddt_array
        self.min_plddt = min_plddt
        self._residue_counter = 0  # tracks 0-indexed residue position

    def accept_residue(self, residue):
        idx = self._residue_counter
        self._residue_counter += 1
        if idx < len(self.plddt):
            return float(self.plddt[idx]) >= self.min_plddt
        return True  # Accept if out of range (safety fallback)


def prepare_predicted_structure(
    pdb_path: str,
    plddt_array: np.ndarray,
    output_path: str,
    min_plddt: float = 70.0,
    remove_disordered_loops: bool = True,
) -> dict:
    """
    Prepare an AlphaFold predicted structure for molecular docking.

    Steps:
    1. Filter out low-confidence residues (pLDDT < min_plddt)
    2. Remove HETATM and water records
    3. Save cleaned PDB for downstream processing with PDBFixer/OpenMM

    Args:
        pdb_path: Input AlphaFold PDB path
        plddt_array: Per-residue pLDDT values
        output_path: Output cleaned PDB path
        min_plddt: Minimum pLDDT for residue retention
        remove_disordered_loops: If True, remove residues with pLDDT < 50

    Returns:
        Dict with statistics about removed residues.
    """
    parser = PDBParser(QUIET=True)
    structure = parser.get_structure("af", pdb_path)

    n_total = sum(1 for _ in structure[0].get_residues())
    threshold = min_plddt if remove_disordered_loops else 50.0

    selector = HighConfidenceSelect(plddt_array, min_plddt=threshold)

    io = PDBIO()
    io.set_structure(structure)
    io.save(output_path, selector)

    n_kept = selector._residue_counter  # Approximation
    n_removed = n_total - n_kept

    return {
        "input_residues": n_total,
        "residues_kept": n_kept,
        "residues_removed_low_confidence": n_removed,
        "output_path": output_path,
    }

AlphaFold vs ESMFold vs Crystal Structure: When to Use What

Pocket Detection on Predicted Structures

import subprocess
import json
from pathlib import Path

def run_fpocket(pdb_path: str, output_dir: str) -> list[dict]:
    """
    Run fpocket on a structure and parse the top pockets.

    Args:
        pdb_path: Path to input PDB (AlphaFold or experimental)
        output_dir: Directory for fpocket output

    Returns:
        List of pocket dicts sorted by druggability score (descending).
    """
    Path(output_dir).mkdir(parents=True, exist_ok=True)
    result = subprocess.run(
        ["fpocket", "-f", pdb_path, "-o", output_dir],
        capture_output=True, text=True,
    )

    pockets = []
    info_file = Path(output_dir) / (Path(pdb_path).stem + "_info.txt")
    if not info_file.exists():
        return []

    current_pocket = {}
    with open(info_file) as f:
        for line in f:
            line = line.strip()
            if line.startswith("Pocket"):
                if current_pocket:
                    pockets.append(current_pocket)
                pocket_num = int(line.split()[1])
                current_pocket = {"pocket_id": pocket_num}
            elif ":" in line:
                key, val = line.split(":", 1)
                try:
                    current_pocket[key.strip().lower().replace(" ", "_")] = float(val.strip())
                except ValueError:
                    current_pocket[key.strip().lower().replace(" ", "_")] = val.strip()
    if current_pocket:
        pockets.append(current_pocket)

    # Sort by druggability score
    pockets.sort(
        key=lambda p: p.get("druggability_score", 0.0),
        reverse=True,
    )
    return pockets


def filter_pockets_by_plddt(
    pockets: list[dict],
    pocket_residue_plddt: dict[int, float],
    min_plddt: float = 70.0,
) -> list[dict]:
    """
    Filter fpocket hits to those whose lining residues are high-confidence.

    Args:
        pockets: List of fpocket pocket dicts
        pocket_residue_plddt: {residue_idx: plddt_value} for each pocket-lining residue
        min_plddt: Minimum mean pLDDT for pocket acceptance

    Returns:
        Filtered list of reliable pockets.
    """
    reliable = []
    for pocket in pockets:
        residue_ids = pocket.get("residue_ids", [])
        if not residue_ids:
            reliable.append(pocket)  # Accept if no residue info available
            continue
        plddt_vals = [pocket_residue_plddt.get(r, 100.0) for r in residue_ids]
        mean_plddt = sum(plddt_vals) / len(plddt_vals)
        if mean_plddt >= min_plddt:
            pocket["mean_pocket_plddt"] = round(mean_plddt, 1)
            reliable.append(pocket)
    return reliable

Key Takeaways

• Always gate AlphaFold structures by pLDDT before docking: require mean pLDDT >70 for the pocket-lining residues. Docking into low-pLDDT regions produces meaningless results.
• For known targets, AlphaFold2 structures are often superior to apo crystal structures for virtual screening because AlphaFold captures the sequence-encoded preferred conformation, which may be closer to the ligand-bound form.
• PAE matrices from AlphaFold3/Multimer are essential for PPI and ternary complex modelling — low inter-chain PAE confirms reliable domain-domain orientations for docking.
• ESMFold (single-sequence, ~1 second per protein) is appropriate for large-scale variant scanning. For any primary docking campaign, use AlphaFold2/3 with MSA for higher accuracy.
• For targets with disordered activation loops (kinases, GPCRs), predicted structures should be complemented with MD sampling or ensemble docking to capture multiple pocket conformations.

Try it live

Explore a related model hands-on: open Protein Folding Visualiser and change its parameters while it runs. Everything stays in your browser.