蛋白质结构预测实战（AlphaFold2 / ESMFold / ColabFold）¶

# FASTA格式示例
cat > target_protein.fasta << 'EOF'
>protein_target
MVLSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYFPHFDLSH
GSAQVKGHGKKVADALTNAVAHVDDMPNALSALSDLHAHKLRVDPVNFKLL
SHCLLVTLAAHLPAEFTPAVHASLDKFLASVSTVLTSKYR
EOF

# 从UniProt下载序列
curl -s "https://rest.uniprot.org/uniprotkb/P69905.fasta" > P69905.fasta

# 检查序列长度
grep -v "^>" target_protein.fasta | tr -d '\n' | wc -c

# AlphaFold2默认MSA流程（内部自动调用）
# 数据库路径
DB_DIR="/data/alphafold_databases"

# JackHMMER搜索UniRef90
jackhmmer \
  --cpu 8 \
  -N 1 \
  --noali \
  --F1 0.0005 --F2 0.00005 --F3 0.0000005 \
  --incdomE 0.0001 \
  -A output.sto \
  target_protein.fasta \
  ${DB_DIR}/uniref90/uniref90.fasta

# ColabFold使用MMseqs2（快速模式）
colabfold_search \
  --db-load-mode 2 \
  --threads 16 \
  target_protein.fasta \
  ${DB_DIR}/colabfold_envdb \
  msa_results/

# HHsearch模板搜索（AlphaFold2内部流程）
hhsearch \
  -i target.a3m \
  -d ${DB_DIR}/pdb70/pdb70 \
  -o template_hits.hhr \
  -maxseq 1000000 \
  -p 20 \
  -Z 500 \
  -z 1 \
  -B 500 \
  -b 1 \
  -loc \
  -mact 0.35

# AlphaFold2本地运行（Docker方式）
python3 docker/run_docker.py \
  --fasta_paths=target_protein.fasta \
  --max_template_date=2024-01-01 \
  --model_preset=monomer_ptm \
  --db_preset=full_dbs \
  --data_dir=/data/alphafold_databases \
  --output_dir=af2_output/ \
  --gpu_devices=0

# 非Docker方式（直接运行）
python3 run_alphafold.py \
  --fasta_paths=target_protein.fasta \
  --output_dir=af2_output/ \
  --data_dir=/data/alphafold_databases \
  --uniref90_database_path=${DB_DIR}/uniref90/uniref90.fasta \
  --mgnify_database_path=${DB_DIR}/mgnify/mgy_clusters_2022_05.fa \
  --template_mmcif_dir=${DB_DIR}/pdb_mmcif/mmcif_files \
  --obsolete_pdbs_path=${DB_DIR}/pdb_mmcif/obsolete.dat \
  --bfd_database_path=${DB_DIR}/bfd/bfd_metaclust_clu_complete_id30_c90_final_seq.sorted_opt \
  --uniref30_database_path=${DB_DIR}/uniref30/UniRef30_2021_03 \
  --pdb70_database_path=${DB_DIR}/pdb70/pdb70 \
  --max_template_date=2024-01-01 \
  --model_preset=monomer_ptm \
  --db_preset=full_dbs \
  --use_gpu_relax=True

# ESMFold Python API
import torch
import esm

# 加载模型
model = esm.pretrained.esmfold_v1()
model = model.eval().cuda()

# 可选：减少内存使用
model.set_chunk_size(128)

# 预测
sequence = "MVLSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYFPHFDLSH"

with torch.no_grad():
    output = model.infer_pdb(sequence)

# 保存PDB文件
with open("esmfold_prediction.pdb", "w") as f:
    f.write(output)

# 获取pLDDT分数
with torch.no_grad():
    output = model.infer(sequence)
    plddt = output["plddt"][0, :, 1]  # per-residue pLDDT
    mean_plddt = plddt.mean().item()
    print(f"Mean pLDDT: {mean_plddt:.2f}")

# 使用ESMFold命令行（通过esm包）
python -c "
import esm
import torch

model = esm.pretrained.esmfold_v1()
model = model.eval().cuda()

with open('target_protein.fasta') as f:
    lines = f.readlines()
    seq = ''.join(l.strip() for l in lines if not l.startswith('>'))

with torch.no_grad():
    pdb_str = model.infer_pdb(seq)

with open('esmfold_output.pdb', 'w') as f:
    f.write(pdb_str)
print('ESMFold prediction saved.')
"

# LocalColabFold安装
# 前提：已安装conda
wget https://raw.githubusercontent.com/YoshitakaMo/localcolabfold/main/install_colabfold_linux.sh
bash install_colabfold_linux.sh

# 将colabfold加入PATH
export PATH="/path/to/localcolabfold/colabfold-conda/bin:$PATH"

# 运行ColabFold批量预测
colabfold_batch \
  target_protein.fasta \
  colabfold_output/ \
  --num-recycle 3 \
  --num-models 5 \
  --amber \
  --templates \
  --num-seeds 1 \
  --use-gpu-relax

# 单序列模式（类似ESMFold场景）
colabfold_batch \
  target_protein.fasta \
  colabfold_output_single/ \
  --msa-mode single_sequence \
  --num-recycle 6

# 使用自定义MSA
colabfold_batch \
  custom_msa.a3m \
  colabfold_output_custom/ \
  --num-recycle 3 \
  --num-models 5

import json
import numpy as np
import matplotlib.pyplot as plt

# 读取AlphaFold2输出的置信度文件
with open("af2_output/ranking_debug.json") as f:
    ranking = json.load(f)
print("模型排名:", ranking["order"])
print("pLDDT分数:", {k: f"{v:.1f}" for k,v in ranking["plddts"].items()})

# 读取PAE矩阵
with open("af2_output/result_model_1_ptm.pkl", "rb") as f:
    import pickle
    result = pickle.load(f)

pae = result["predicted_aligned_error"]
plddt = result["plddt"]

# 绘制pLDDT per-residue图
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# pLDDT图
axes[0].plot(range(1, len(plddt)+1), plddt, color="steelblue", linewidth=0.8)
axes[0].axhline(y=70, color="orange", linestyle="--", label="pLDDT=70")
axes[0].axhline(y=90, color="green", linestyle="--", label="pLDDT=90")
axes[0].set_xlabel("残基位置")
axes[0].set_ylabel("pLDDT")
axes[0].set_title("Per-residue Confidence (pLDDT)")
axes[0].set_ylim(0, 100)
axes[0].legend()

# PAE图
im = axes[1].imshow(pae, cmap="Greens_r", vmin=0, vmax=30)
axes[1].set_xlabel("Scored residue")
axes[1].set_ylabel("Aligned residue")
axes[1].set_title("Predicted Aligned Error (PAE)")
plt.colorbar(im, ax=axes[1], label="Expected position error (Å)")

plt.tight_layout()
plt.savefig("quality_assessment.png", dpi=300, bbox_inches="tight")

# 使用OpenMM进行能量最小化
from openmm.app import *
from openmm import *
from openmm.unit import *

# 加载PDB
pdb = PDBFile("af2_prediction_unrelaxed.pdb")
forcefield = ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')

# 创建系统
modeller = Modeller(pdb.topology, pdb.positions)
modeller.addHydrogens(forcefield)
system = forcefield.createSystem(modeller.topology,
    nonbondedMethod=NoCutoff,
    constraints=HBonds)

# 能量最小化
integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
simulation = Simulation(modeller.topology, system, integrator)
simulation.context.setPositions(modeller.positions)

# 运行最小化
simulation.minimizeEnergy(maxIterations=1000)

# 保存
state = simulation.context.getState(getPositions=True)
with open("af2_prediction_relaxed.pdb", "w") as f:
    PDBFile.writeFile(simulation.topology, state.getPositions(), f)

# PyMOL脚本：按pLDDT着色
cat > color_by_plddt.pml << 'EOF'
load af2_prediction.pdb, prediction

# 按B-factor（存储pLDDT）着色
spectrum b, blue_white_red, prediction, minimum=50, maximum=100

# 设置显示
bg_color white
set cartoon_fancy_helices, 1
set cartoon_side_chain_helper, 1
set ray_opaque_background, 1

# 保存图像
ray 2400, 1800
png prediction_colored.png, dpi=300
EOF

pymol -c color_by_plddt.pml

# 使用TM-align比较两个结构
TMalign af2_prediction.pdb experimental_structure.pdb -o alignment.sup

# AlphaFold-Multimer输入格式（两条链）
cat > complex.fasta << 'EOF'
>chain_A
MVLSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYFPHFDLSH
>chain_B  
MHSSIVLATVLFVAIASASKTRELCMKSLEHAKVGTSKEAKQDGIDLYKHMFE
EOF

# 运行AlphaFold-Multimer
python3 run_alphafold.py \
  --fasta_paths=complex.fasta \
  --output_dir=multimer_output/ \
  --data_dir=/data/alphafold_databases \
  --model_preset=multimer \
  --max_template_date=2024-01-01 \
  --use_gpu_relax=True

# ColabFold多聚体预测
colabfold_batch \
  complex.fasta \
  colabfold_multimer_output/ \
  --num-recycle 6 \
  --num-models 5 \
  --amber

#!/bin/bash
# ==============================================
# 蛋白质结构预测完整流程
# 综合使用AlphaFold2、ESMFold和ColabFold
# ==============================================

set -euo pipefail

TARGET="target_protein"
FASTA="${TARGET}.fasta"
OUTDIR="structure_predictions"
AF2_DB="/data/alphafold_databases"

mkdir -p ${OUTDIR}/{alphafold2,esmfold,colabfold,comparison}

echo "=== 目标蛋白信息 ==="
SEQ_LEN=$(grep -v "^>" ${FASTA} | tr -d '\n' | wc -c)
echo "序列长度: ${SEQ_LEN} 残基"

# ---- 方法1: ColabFold（推荐首选） ----
echo ""
echo "=== ColabFold预测 ==="
colabfold_batch \
  ${FASTA} \
  ${OUTDIR}/colabfold/ \
  --num-recycle 3 \
  --num-models 5 \
  --amber \
  --templates \
  --num-seeds 1 \
  --use-gpu-relax

# ---- 方法2: AlphaFold2（最高精度） ----
echo ""
echo "=== AlphaFold2预测 ==="
python3 run_alphafold.py \
  --fasta_paths=${FASTA} \
  --output_dir=${OUTDIR}/alphafold2/ \
  --data_dir=${AF2_DB} \
  --model_preset=monomer_ptm \
  --db_preset=full_dbs \
  --max_template_date=2024-01-01 \
  --use_gpu_relax=True

# ---- 方法3: ESMFold（快速） ----
echo ""
echo "=== ESMFold预测 ==="
python3 << 'PYTHON'
import torch
import esm

model = esm.pretrained.esmfold_v1()
model = model.eval().cuda()
model.set_chunk_size(128)

with open("${FASTA}") as f:
    lines = f.readlines()
    seq = ''.join(l.strip() for l in lines if not l.startswith('>'))

print(f"预测序列长度: {len(seq)}")

with torch.no_grad():
    pdb_str = model.infer_pdb(seq)

with open("${OUTDIR}/esmfold/prediction.pdb", "w") as f:
    f.write(pdb_str)

# 获取pLDDT
with torch.no_grad():
    output = model.infer(seq)
    plddt = output["plddt"][0, :, 1].cpu().numpy()
    print(f"Mean pLDDT: {plddt.mean():.2f}")
    print(f"Median pLDDT: {float(sorted(plddt)[len(plddt)//2]):.2f}")

PYTHON

echo ""
echo "=== 结构比较 ==="
# 使用TM-align比较三种方法的结果
TMalign \
  ${OUTDIR}/colabfold/*_relaxed_rank_001*.pdb \
  ${OUTDIR}/alphafold2/${TARGET}/ranked_0.pdb \
  -o ${OUTDIR}/comparison/cf_vs_af2

TMalign \
  ${OUTDIR}/esmfold/prediction.pdb \
  ${OUTDIR}/alphafold2/${TARGET}/ranked_0.pdb \
  -o ${OUTDIR}/comparison/esm_vs_af2

echo "=== 完成 ==="

#!/usr/bin/env python3
"""
蛋白质结构预测质量评估和可视化脚本
"""

import json
import pickle
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path
from Bio.PDB import PDBParser, DSSP, Superimposer

def load_af2_results(output_dir):
    """加载AlphaFold2预测结果"""
    results = {}
    output_path = Path(output_dir)

    # 读取排名
    with open(output_path / "ranking_debug.json") as f:
        ranking = json.load(f)
    results["ranking"] = ranking

    # 读取最佳模型的详细结果
    best_model = ranking["order"][0]
    pkl_file = output_path / f"result_{best_model}.pkl"
    with open(pkl_file, "rb") as f:
        model_result = pickle.load(f)

    results["plddt"] = model_result["plddt"]
    results["pae"] = model_result.get("predicted_aligned_error", None)
    results["ptm"] = model_result.get("ptm", None)

    return results

def plot_quality(results, output_prefix):
    """绘制质量评估图"""
    fig, axes = plt.subplots(1, 3, figsize=(18, 5))

    plddt = results["plddt"]

    # 1. pLDDT per-residue
    residues = np.arange(1, len(plddt)+1)
    axes[0].plot(residues, plddt, linewidth=0.8, color="steelblue")
    axes[0].fill_between(residues, 0, plddt, alpha=0.3, color="steelblue")
    axes[0].axhline(70, color="orange", linestyle="--", alpha=0.7, label="pLDDT=70")
    axes[0].axhline(90, color="green", linestyle="--", alpha=0.7, label="pLDDT=90")
    axes[0].set_xlabel("Residue Index")
    axes[0].set_ylabel("pLDDT")
    axes[0].set_title(f"Per-residue Confidence\n(Mean: {np.mean(plddt):.1f})")
    axes[0].set_ylim(0, 100)
    axes[0].legend()

    # 2. pLDDT分布直方图
    axes[1].hist(plddt, bins=50, color="steelblue", edgecolor="white")
    axes[1].axvline(70, color="orange", linestyle="--")
    axes[1].axvline(90, color="green", linestyle="--")
    axes[1].set_xlabel("pLDDT")
    axes[1].set_ylabel("Count")
    axes[1].set_title("pLDDT Distribution")

    # 3. PAE热图
    if results["pae"] is not None:
        im = axes[2].imshow(results["pae"], cmap="Greens_r", vmin=0, vmax=30)
        axes[2].set_xlabel("Scored residue")
        axes[2].set_ylabel("Aligned residue")
        axes[2].set_title("Predicted Aligned Error")
        plt.colorbar(im, ax=axes[2], label="Error (Å)")

    plt.tight_layout()
    plt.savefig(f"{output_prefix}_quality.png", dpi=300, bbox_inches="tight")
    plt.close()

    # 打印汇总
    print(f"Mean pLDDT: {np.mean(plddt):.1f}")
    print(f"Residues with pLDDT > 90: {np.sum(plddt > 90)}/{len(plddt)} ({100*np.mean(plddt>90):.1f}%)")
    print(f"Residues with pLDDT > 70: {np.sum(plddt > 70)}/{len(plddt)} ({100*np.mean(plddt>70):.1f}%)")
    print(f"Residues with pLDDT < 50: {np.sum(plddt < 50)}/{len(plddt)} ({100*np.mean(plddt<50):.1f}%)")
    if results["ptm"] is not None:
        print(f"pTM score: {results['ptm']:.3f}")

if __name__ == "__main__":
    results = load_af2_results("structure_predictions/alphafold2/target_protein/")
    plot_quality(results, "structure_predictions/comparison/af2")