药物靶点预测分析（分子对接/药效团/ADMET预测/网络药理学）¶

# ===== 从PDB下载蛋白结构 =====
# 下载PDB文件
wget https://files.rcsb.org/download/6LU7.pdb
# 6LU7是SARS-CoV-2 Main protease (Mpro)

# ===== AlphaFold预测结构 =====
# 如果PDB中没有目标蛋白结构
# 从AlphaFold DB下载
wget https://alphafold.ebi.ac.uk/files/AF-P0DTD1-F1-model_v4.pdb

# ===== 使用PyMOL准备蛋白 =====
# pymol命令
# load 6LU7.pdb
# remove solvent          # 去水
# remove resn NAG+FUC+MAN  # 去糖基
# h_add                    # 加氢
# save receptor_prepared.pdb

# ===== 使用Open Babel准备 =====
# 安装
conda install -c conda-forge openbabel

# PDB转PDBQT（AutoDock Vina需要的格式）
# 使用MGLTools/AutoDockTools更专业
# 安装ADFR suite (含prepare_receptor)
wget https://ccsb.scripps.edu/adfr/download/1038/
tar -xzf ADFRsuite_x86_64Linux_1.0.tar.gz

# 准备受体（去水、加氢、加电荷、转PDBQT）
prepare_receptor -r 6LU7.pdb -o receptor.pdbqt -A hydrogens

# ===== Python准备蛋白（使用MDAnalysis/BioPython） =====
from Bio.PDB import PDBParser, PDBIO, Select
import subprocess

# 解析PDB
parser = PDBParser()
structure = parser.get_structure("receptor", "6LU7.pdb")

# 只保留蛋白链（去除水和配体）
class ProteinOnly(Select):
    def accept_residue(self, residue):
        return residue.id[0] == " "  # 只保留标准残基

io = PDBIO()
io.set_structure(structure)
io.save("receptor_clean.pdb", ProteinOnly())

# 使用Meeko准备受体（推荐用于Vina）
# pip install meeko
from meeko import MoleculePreparation, PDBQTWriterLegacy

# 或使用命令行
subprocess.run([
    "prepare_receptor",
    "-r", "receptor_clean.pdb",
    "-o", "receptor.pdbqt",
    "-A", "hydrogens"
])

# ===== 配体准备 =====
# 从PubChem下载SDF格式
wget "https://pubchem.ncbi.nlm.nih.gov/rest/pug/compound/name/quercetin/SDF" -O quercetin.sdf

# SDF转PDBQT（使用Open Babel）
obabel quercetin.sdf -O quercetin.pdbqt --gen3d -h

# 或使用Meeko（推荐，更好的扭转角处理）
mk_prepare_ligand.py -i quercetin.sdf -o quercetin.pdbqt

# ===== AutoDock Vina对接 =====
# 安装Vina
conda install -c conda-forge vina
# 或下载二进制文件

# 配置文件 (config.txt)
cat > config.txt << 'EOF'
receptor = receptor.pdbqt
ligand = quercetin.pdbqt
center_x = -10.0
center_y = 15.0
center_z = 68.0
size_x = 20
size_y = 20
size_z = 20
exhaustiveness = 32
num_modes = 10
energy_range = 3
EOF

# 运行对接
vina --config config.txt --out docked_output.pdbqt --log docking_log.txt

# 查看结果
cat docking_log.txt
# mode |   affinity | dist from best mode
#      | (kcal/mol) | rmsd l.b.| rmsd u.b.
# -----+------------+----------+----------
#    1       -7.8         0.000      0.000
#    2       -7.5         1.234      2.456
#    3       -7.2         2.345      4.567

# ===== Python Vina接口（推荐） =====
from vina import Vina

# 创建Vina对象
v = Vina(sf_name="vina")

# 设置受体
v.set_receptor("receptor.pdbqt")

# 设置配体
v.set_ligand_from_file("quercetin.pdbqt")

# 设置对接盒子（基于共结晶配体位置或活性位点中心）
v.compute_vina_maps(
    center=[-10.0, 15.0, 68.0],
    box_size=[20, 20, 20]
)

# 对接
v.dock(exhaustiveness=32, n_poses=10)

# 获取结果
energies = v.energies()
print(f"Best binding energy: {energies[0][0]:.2f} kcal/mol")

# 保存对接构象
v.write_poses("docked_poses.pdbqt", n_poses=5, overwrite=True)

# ===== 批量对接（虚拟筛选） =====
import os
import pandas as pd

ligand_dir = "ligands_pdbqt/"
results = []

for lig_file in os.listdir(ligand_dir):
    if lig_file.endswith(".pdbqt"):
        v = Vina(sf_name="vina")
        v.set_receptor("receptor.pdbqt")
        v.set_ligand_from_file(os.path.join(ligand_dir, lig_file))
        v.compute_vina_maps(center=[-10, 15, 68], box_size=[20, 20, 20])
        v.dock(exhaustiveness=16, n_poses=5)

        energy = v.energies()[0][0]
        results.append({
            "ligand": lig_file.replace(".pdbqt", ""),
            "binding_energy": energy
        })

# 排序结果
df = pd.DataFrame(results).sort_values("binding_energy")
df.to_csv("virtual_screening_results.csv", index=False)
print(df.head(20))

# ===== PLIP分析互作 =====
# Protein-Ligand Interaction Profiler
# pip install plip

from plip.structure.preparation import PDBComplex

# 准备复合物文件（蛋白+对接配体）
# 先用PyMOL或obabel合并蛋白和对接pose
my_mol = PDBComplex()
my_mol.load_pdb("complex.pdb")
my_mol.analyze()

for bsid, site in sorted(my_mol.interaction_sets.items()):
    print(f"\n=== Binding Site: {bsid} ===")
    print(f"Hydrogen bonds: {len(site.hbonds_ldon + site.hbonds_pdon)}")
    print(f"Hydrophobic contacts: {len(site.hydrophobic_contacts)}")
    print(f"Pi-stacking: {len(site.pistacking)}")
    print(f"Salt bridges: {len(site.saltbridge_lneg + site.saltbridge_pneg)}")

    for hbond in site.hbonds_pdon + site.hbonds_ldon:
        print(f"  H-bond: {hbond.resnr}{hbond.restype} - distance: {hbond.distance_ad:.2f}Å")

# ===== PyMOL可视化脚本 =====
pymol_script = """
load receptor.pdb, protein
load docked_poses.pdbqt, ligand
split_states ligand
delete ligand
# 只保留最佳构象
hide everything
show cartoon, protein
show sticks, ligand_0001
color cyan, ligand_0001

# 显示结合位点残基
select binding_site, protein within 4 of ligand_0001
show sticks, binding_site
color yellow, binding_site

# 显示氢键
distance hbonds, ligand_0001, binding_site, 3.5, mode=2

# 表面
show surface, protein
set surface_color, white
set transparency, 0.7

# 保存图片
ray 2400, 2400
png docking_result.png, dpi=300
"""

with open("visualize_docking.pml", "w") as f:
    f.write(pymol_script)
# pymol visualize_docking.pml

# ===== SwissADME (在线工具API化) =====
# http://www.swissadme.ch/
# 输入SMILES，获取ADME参数

# ===== ADMETlab 2.0 =====
# https://admetmesh.scbdd.com/

# ===== RDKit计算基本描述符 =====
from rdkit import Chem
from rdkit.Chem import Descriptors, Lipinski, Draw
import pandas as pd

# 计算Lipinski类药五规则
def lipinski_filter(smiles):
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return None

    mw = Descriptors.ExactMolWt(mol)
    logp = Descriptors.MolLogP(mol)
    hbd = Descriptors.NumHDonors(mol)
    hba = Descriptors.NumHAcceptors(mol)
    tpsa = Descriptors.TPSA(mol)
    rotatable = Descriptors.NumRotatableBonds(mol)

    violations = 0
    if mw > 500: violations += 1
    if logp > 5: violations += 1
    if hbd > 5: violations += 1
    if hba > 10: violations += 1

    return {
        "MW": mw,
        "LogP": logp,
        "HBD": hbd,
        "HBA": hba,
        "TPSA": tpsa,
        "RotatableBonds": rotatable,
        "Lipinski_violations": violations,
        "Drug_like": violations <= 1
    }

# 批量计算
compounds = {
    "Quercetin": "O=C1C(O)=C(OC2=CC(O)=CC(O)=C12)C1=CC(O)=C(O)C=C1",
    "Aspirin": "CC(=O)OC1=CC=CC=C1C(O)=O",
    "Metformin": "CN(C)C(=N)NC(=N)N"
}

results = []
for name, smiles in compounds.items():
    props = lipinski_filter(smiles)
    if props:
        props["Name"] = name
        props["SMILES"] = smiles
        results.append(props)

df_admet = pd.DataFrame(results)
print(df_admet[["Name", "MW", "LogP", "HBD", "HBA", "TPSA", "Drug_like"]])

# ===== pkCSM在线预测 =====
# http://biosig.unimelb.edu.au/pkcsm/
# 支持批量SMILES输入，预测完整ADMET参数

# ===== DeepPurpose（深度学习ADMET） =====
# pip install DeepPurpose
# 使用预训练模型预测药物-靶点相互作用

# ===== 网络药理学R语言流程 =====
library(clusterProfiler)
library(org.Hs.eg.db)
library(ggplot2)

# 步骤1：获取药物靶点
# 来源：SwissTargetPrediction, PharmMapper, SEA, STITCH
# 假设已获取药物靶点列表
drug_targets <- read.csv("drug_targets.csv")  # 列：compound, target

# 步骤2：获取疾病靶点
# 来源：GeneCards, OMIM, DisGeNET, TTD
disease_targets <- read.csv("disease_targets.csv")  # 列：gene, score

# 步骤3：取交集靶点
common_targets <- intersect(drug_targets$target, disease_targets$gene)
cat("共同靶点数:", length(common_targets), "\n")

# 步骤4：PPI网络构建（参见蛋白互作网络分析）
library(STRINGdb)
string_db <- STRINGdb$new(version = "12.0", species = 9606, score_threshold = 400)
mapped <- string_db$map(data.frame(gene = common_targets), "gene", removeUnmappedRows = TRUE)
ppi <- string_db$get_interactions(mapped$STRING_id)

# 步骤5：GO和KEGG富集分析
ego <- enrichGO(
  gene = common_targets,
  OrgDb = org.Hs.eg.db,
  keyType = "SYMBOL",
  ont = "BP",
  pAdjustMethod = "BH",
  qvalueCutoff = 0.05
)

ekegg <- enrichKEGG(
  gene = bitr(common_targets, fromType = "SYMBOL", toType = "ENTREZID",
              OrgDb = org.Hs.eg.db)$ENTREZID,
  organism = "hsa",
  pAdjustMethod = "BH"
)

# 可视化
dotplot(ego, showCategory = 20, title = "GO Enrichment")
dotplot(ekegg, showCategory = 20, title = "KEGG Enrichment")

# 步骤6：构建药物-靶点-通路网络
# 导出到Cytoscape
network_edges <- data.frame(
  source = c(rep("Drug_A", length(common_targets)),
             rep(common_targets, each = 1)),
  target = c(common_targets,
             rep("Disease_X", length(common_targets))),
  type = c(rep("drug-target", length(common_targets)),
           rep("target-disease", length(common_targets)))
)
write.csv(network_edges, "network_for_cytoscape.csv", row.names = FALSE)

# ===== RDKit药效团生成 =====
from rdkit import Chem
from rdkit.Chem import AllChem, Draw
from rdkit.Chem.Pharm3D import Pharmacophore, EmbedLib
from rdkit.Chem import ChemicalFeatures
from rdkit import RDConfig
import os

# 加载药效团特征定义
fdefName = os.path.join(RDConfig.RDDataDir, 'BaseFeatures.fdef')
featFactory = ChemicalFeatures.BuildFeatureFactory(fdefName)

# 读取活性分子
mol = Chem.MolFromSmiles("O=C1C(O)=C(OC2=CC(O)=CC(O)=C12)C1=CC(O)=C(O)C=C1")
mol = Chem.AddHs(mol)
AllChem.EmbedMolecule(mol, AllChem.ETKDGv3())

# 提取药效团特征
feats = featFactory.GetFeaturesForMol(mol)
for feat in feats:
    print(f"{feat.GetFamily()} at position {feat.GetPos()}")
    # 输出如：Donor, Acceptor, Hydrophobe, Aromatic, etc.

# ===== 多分子药效团对齐 =====
# 需要多个已知活性分子来建立共同药效团
from rdkit.Chem import rdMolAlign

active_smiles = [
    "O=C1C(O)=C(OC2=CC(O)=CC(O)=C12)C1=CC(O)=C(O)C=C1",  # Quercetin
    "O=C1C(O)=C(OC2=CC(O)=CC(O)=C12)C1=CC=C(O)C=C1",      # Kaempferol
    "OC1=CC(O)=C2C(=O)C(O)=C(OC2=C1)C1=CC(O)=C(O)C(O)=C1" # Myricetin
]

mols = [Chem.MolFromSmiles(s) for s in active_smiles]
mols = [Chem.AddHs(m) for m in mols]
[AllChem.EmbedMolecule(m, AllChem.ETKDGv3()) for m in mols]

# 对齐分子
ref_mol = mols[0]
for mol in mols[1:]:
    AllChem.AlignMol(mol, ref_mol)

#!/bin/bash
# ============================================
# 药物靶点预测完整流程
# ============================================

# 1. 安装必要工具
conda create -n drug_discovery python=3.10
conda activate drug_discovery
pip install vina meeko rdkit numpy pandas
conda install -c conda-forge openbabel pymol-open-source

# 2. 下载蛋白结构
wget https://files.rcsb.org/download/6LU7.pdb -O receptor.pdb

# 3. 准备受体
prepare_receptor -r receptor.pdb -o receptor.pdbqt -A hydrogens

# 4. 准备配体（批量）
mkdir -p ligands_pdbqt
for sdf in ligands_sdf/*.sdf; do
  name=$(basename $sdf .sdf)
  mk_prepare_ligand.py -i $sdf -o ligands_pdbqt/${name}.pdbqt
done

# 5. 批量对接
mkdir -p docking_results
for lig in ligands_pdbqt/*.pdbqt; do
  name=$(basename $lig .pdbqt)
  vina --receptor receptor.pdbqt \
       --ligand $lig \
       --center_x -10 --center_y 15 --center_z 68 \
       --size_x 20 --size_y 20 --size_z 20 \
       --exhaustiveness 32 \
       --out docking_results/${name}_docked.pdbqt \
       --log docking_results/${name}_log.txt
done

# 6. 汇总结果
echo "Ligand,Binding_Energy" > docking_summary.csv
for log in docking_results/*_log.txt; do
  name=$(basename $log _log.txt)
  energy=$(grep "^   1" $log | awk '{print $2}')
  echo "$name,$energy" >> docking_summary.csv
done
sort -t',' -k2 -n docking_summary.csv