154_蛋白质互作预测¶

# === AlphaFold-Multimer 本地安装 ===
# 需要: NVIDIA GPU (A100推荐), ~2TB磁盘(数据库), Docker

# 1. 克隆AlphaFold仓库
git clone https://github.com/google-deepmind/alphafold.git
cd alphafold

# 2. 下载数据库 (约2.6TB)
bash scripts/download_all_data.sh /path/to/databases/

# 3. 准备输入FASTA文件 (多链用">"分隔)
cat > complex.fasta << 'EOF'
>ChainA
MVLSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYFPHFDLSH
GSAQVKGHGKKVADALTNAVAHVDDMPNALSALSDLHAHKLRVDPVNFKLL
>ChainB
MVHLTPEEKSAVTALWGKVNVDEVGGEALGRLLVVYPWTQRFFESFGDLST
PDAVMGNPKVKAHGKKVLGAFSDGLAHLDNLKGTFATLSELHCDKLHVDPE
EOF

# 4. 运行AlphaFold-Multimer
python run_alphafold.py \
  --fasta_paths=complex.fasta \
  --model_preset=multimer \
  --data_dir=/path/to/databases/ \
  --output_dir=output/ \
  --max_template_date=2024-01-01 \
  --db_preset=full_dbs \
  --use_gpu_relax=True

# 输出文件:
# output/complex/
#   ranked_0.pdb          : 最佳模型
#   ranked_1.pdb          : 次佳模型
#   ...
#   ranking_debug.json    : 排名分数
#   result_model_*.pkl    : 包含pLDDT, PAE等信息

# === 安装ColabFold本地版 ===
pip install colabfold[alphafold]
# 或使用conda
conda install -c conda-forge -c bioconda colabfold

# === 运行ColabFold预测复合物 ===
# 多链序列用":"分隔写在一个条目中
cat > complex.fasta << 'EOF'
>complex_AB
MVLSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYFPHFDLSH:MVHLTPEEKSAVTALWGKVNVDEVGGEALGRLLVVYPWTQRFFESFGDLST
EOF

# 使用colabfold_batch
colabfold_batch \
  complex.fasta \
  output_dir/ \
  --num-recycle 3 \
  --num-models 5 \
  --amber \
  --use-gpu-relax \
  --model-type alphafold2_multimer_v3

# 参数说明:
# --num-recycle: recycling次数(提高准确度, 默认3)
# --num-models: 使用的模型数(1-5)
# --amber: 是否做AMBER力场relaxation
# --model-type: 模型类型(multimer_v3最新)

# === 结果文件 ===
# output_dir/
#   complex_AB_unrelaxed_rank_001.pdb  : 最佳未松弛模型
#   complex_AB_relaxed_rank_001.pdb    : AMBER松弛后模型
#   complex_AB_scores_rank_001.json    : 分数(pLDDT, PAE, iPTM)
#   complex_AB_coverage.png            : MSA覆盖度图
#   complex_AB_pae.png                 : PAE热图
#   complex_AB_plddt.png               : pLDDT分布图

# === 使用MMseqs2 API (最快方式) ===
# ColabFold默认使用MMseqs2服务器进行MSA搜索
# 如需本地MSA搜索:
colabfold_search \
  complex.fasta \
  /path/to/colabfold_db/ \
  msas/ \
  --threads 16

colabfold_batch \
  msas/ \
  output_dir/ \
  --num-models 5

# === HDOCK 网页版使用 ===
# 网址: http://hdock.phys.hust.edu.cn/

# === HDOCK 本地版 ===
# 下载: http://hdock.phys.hust.edu.cn/software.php
# 需要申请license

# 1. 准备输入文件
# receptor.pdb : 受体蛋白结构
# ligand.pdb   : 配体蛋白结构

# 2. 运行HDOCK
hdock receptor.pdb ligand.pdb \
  -out hdock_output \
  -nout 100 \
  -angle 15

# 参数:
# -out: 输出前缀
# -nout: 输出模型数量
# -angle: 旋转角度步长(度)

# 3. 输出
# hdock_output_*.pdb : 对接模型
# Hdock.out : 评分文件

# === 使用约束信息 ===
# 如果知道某些残基在界面上, 可以加约束
# restraint.txt 格式:
# RES A 45 50   (受体A链45-50号残基)
# RES B 100 110 (配体B链100-110号残基)
hdock receptor.pdb ligand.pdb \
  -out hdock_output \
  -rest restraint.txt

import json
import numpy as np
from Bio.PDB import PDBParser, NeighborSearch, Selection
import matplotlib.pyplot as plt

# === 1. AlphaFold质量评估 ===
# 读取scores
with open("output/complex_AB_scores_rank_001.json") as f:
    scores = json.load(f)

# 关键指标
print(f"pTM:  {scores['ptm']:.3f}")    # 单链TM-score预测
print(f"iPTM: {scores['iptm']:.3f}")   # 界面TM-score预测 (>0.8 = 高置信)
print(f"Ranking confidence: {0.8*scores['iptm'] + 0.2*scores['ptm']:.3f}")

# pLDDT分析
plddt = np.array(scores['plddt'])
print(f"Mean pLDDT: {plddt.mean():.1f}")
print(f"pLDDT > 70 fraction: {(plddt > 70).sum() / len(plddt):.2%}")

# PAE热图
pae = np.array(scores['pae'])
plt.figure(figsize=(8, 6))
plt.imshow(pae, cmap='Greens_r', vmin=0, vmax=30)
plt.colorbar(label='PAE (Å)')
plt.xlabel('Residue index (scored)')
plt.ylabel('Residue index (aligned)')
plt.title('Predicted Aligned Error')
plt.savefig('pae_heatmap.png', dpi=150)

# === 2. 界面分析 ===
parser = PDBParser(QUIET=True)
structure = parser.get_structure("complex", "complex_AB_relaxed_rank_001.pdb")
model = structure[0]

# 定义链
chain_A = model['A']
chain_B = model['B']

# 找界面残基 (距离 < 5Å)
atoms_A = Selection.unfold_entities(chain_A, 'A')
atoms_B = Selection.unfold_entities(chain_B, 'A')

ns = NeighborSearch(list(atoms_B))
interface_residues_A = set()
interface_residues_B = set()

for atom in atoms_A:
    close_atoms = ns.search(atom.get_vector().get_array(), 5.0)
    if close_atoms:
        interface_residues_A.add(atom.get_parent().get_id()[1])
        for ca in close_atoms:
            interface_residues_B.add(ca.get_parent().get_id()[1])

print(f"Interface residues in chain A: {len(interface_residues_A)}")
print(f"Interface residues in chain B: {len(interface_residues_B)}")

# === 3. 使用PRODIGY估计结合亲和力 ===
# pip install prodigy
# prodigy complex.pdb --selection A B

# PRODIGY结合亲和力预测
pip install prodigy-prot
prodigy complex_relaxed.pdb --selection A B
# 输出: 预测的Kd和ΔG

# PISA界面分析 (PDBePISA web服务)
# https://www.ebi.ac.uk/pdbe/pisa/

# PyMOL可视化界面
pymol -c << 'EOF'
load complex_relaxed.pdb
select interface_A, chain A and (chain B around 5)
select interface_B, chain B and (chain A around 5)
show sticks, interface_A or interface_B
color red, interface_A
color blue, interface_B
find_pairs interface_A, interface_B, mode=2, cutoff=3.5
save interface_view.png, ray=1
EOF

# === 批量预测蛋白互作 ===
import subprocess
import os
from itertools import combinations

# 蛋白序列列表
proteins = {
    "ProtA": "MVLSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYFPHFDLSH",
    "ProtB": "MVHLTPEEKSAVTALWGKVNVDEVGGEALGRLLVVYPWTQRFFESFGDLST",
    "ProtC": "MNIFEMLRIDEGLRLKIYKDTEGYYTIGIGHLLTKSPSLNAAKSELDKAIGRN",
    "ProtD": "MKWVTFISLLFLFSSAYSRGVFRRDAHKSEVAHRFKDLGEENFKALVLIA"
}

# 生成所有成对组合
pairs = list(combinations(proteins.keys(), 2))

# 为每对生成FASTA并运行ColabFold
for p1, p2 in pairs:
    fasta_content = f">{p1}_{p2}\n{proteins[p1]}:{proteins[p2]}\n"
    fasta_file = f"pairs/{p1}_{p2}.fasta"
    os.makedirs("pairs", exist_ok=True)
    with open(fasta_file, 'w') as f:
        f.write(fasta_content)

# 批量运行
subprocess.run([
    "colabfold_batch",
    "pairs/",
    "results/",
    "--num-recycle", "3",
    "--num-models", "1",  # 快速筛查用1个模型
    "--model-type", "alphafold2_multimer_v3"
])

# === 解析结果, 按iPTM排序 ===
import glob
results = []
for json_file in glob.glob("results/*scores*.json"):
    with open(json_file) as f:
        data = json.load(f)
    pair_name = os.path.basename(json_file).split("_scores")[0]
    results.append({
        "pair": pair_name,
        "iPTM": data["iptm"],
        "pTM": data["ptm"],
        "confidence": 0.8 * data["iptm"] + 0.2 * data["ptm"]
    })

results.sort(key=lambda x: x["confidence"], reverse=True)
for r in results:
    print(f"{r['pair']}: iPTM={r['iPTM']:.3f}, confidence={r['confidence']:.3f}")

# === ColabFold一键安装和运行 ===
pip install colabfold[alphafold]

# 复合物预测
colabfold_batch complex.fasta output/ \
  --num-recycle 3 --num-models 5 --amber

# === 本地AlphaFold Multimer ===
python run_alphafold.py \
  --fasta_paths=complex.fasta \
  --model_preset=multimer \
  --data_dir=/data/af2_db/ \
  --output_dir=af2_output/ \
  --use_gpu_relax

# === HDOCK在线提交 ===
curl -F "receptor=@receptor.pdb" \
     -F "ligand=@ligand.pdb" \
     http://hdock.phys.hust.edu.cn/api/submit

# === PyMOL结构可视化 ===
pymol -c -d "load ranked_0.pdb; spectrum b, blue_white_red; png output.png, ray=1, dpi=300"