AI药物发现与分子生成完全教程
1. AI药物发现概述与技术栈
传统药物研发平均耗时 10-15 年，耗资超过 26 亿美元，且成功率不足 10%。AI 正在从根本上改变这一格局——通过计算方法加速靶点发现、先导化合物筛选和临床试验设计，将部分环节的时间从数年缩短到数周。
AI 药物发现的核心技术栈包括：
- 分子表示学习：将化学分子转化为机器可处理的向量/图结构
- 分子属性预测：预测溶解度、毒性、生物活性等理化性质
- 分子生成与优化：设计具有目标属性的全新分子
- 蛋白质结构预测：预测靶点蛋白的三维结构
- 药物-靶点相互作用：预测分子与蛋白质的结合能力
- 虚拟筛选：从大规模化合物库中高效筛选候选药物
核心工具生态包括 RDKit（化学信息学）、PyTorch Geometric（图神经网络）、DeepChem（药物发现专用框架）和 Hugging Face 的分子模型库。
# 环境准备（核心依赖）
# pip install rdkit-pypi torch-geometric deepchem
from rdkit import Chem
from rdkit.Chem import Draw, Descriptors, AllChem
import numpy as np
# 用RDKit加载一个分子（阿司匹林）
aspirin_smiles = 'CC(=O)Oc1ccccc1C(=O)O'
mol = Chem.MolFromSmiles(aspirin_smiles)
print(f'分子式: {Chem.rdMolDescriptors.CalcMolFormula(mol)}')
print(f'分子量: {Descriptors.MolWt(mol):.2f}')
print(f'LogP: {Descriptors.MolLogP(mol):.2f}')
print(f'氢键供体: {Descriptors.NumHDonors(mol)}')
print(f'氢键受体: {Descriptors.NumHAcceptors(mol)}')
2. 分子表示学习
2.1 SMILES 字符串表示
SMILES（Simplified Molecular Input Line Entry System）是一种将分子结构编码为 ASCII 字符串的线性表示法。它简洁高效，可以作为序列模型（如 Transformer）的输入。
from rdkit import Chem
# 常见分子的SMILES
molecules = {
'阿司匹林': 'CC(=O)Oc1ccccc1C(=O)O',
'布洛芬': 'CC(C)Cc1ccc(cc1)C(C)C(=O)O',
'对乙酰氨基酚': 'CC(=O)Nc1ccc(O)cc1',
'咖啡因': 'Cn1c(=O)c2c(ncn2C)n(C)c1=O',
}
# SMILES验证：解析后重新生成规范化SMILES
for name, smi in molecules.items():
mol = Chem.MolFromSmiles(smi)
canonical = Chem.MolToSmiles(mol)
print(f'{name}: {canonical}')
SELFIES（Self-Referencing Embedded Strings） 是一种更鲁棒的分子表示法，任意字符串都能解码为合法分子，这对生成模型尤为重要。
# pip install selfies
import selfies as sf
# SMILES -> SELFIES
aspirin_smiles = 'CC(=O)Oc1ccccc1C(=O)O'
aspirin_selfies = sf.encoder(aspirin_smiles)
print(f'SELFIES: {aspirin_selfies}')
# SELFIES -> SMILES
decoded = sf.decoder(aspirin_selfies)
print(f'解码回SMILES: {decoded}')
# SELFIES的鲁棒性：随机修改SELFIES仍能解码为合法分子
tokens = list(sf.split_selfies(aspirin_selfies))
tokens[3] = '[C]' # 替换一个token
modified = ''.join(tokens)
decoded_modified = sf.decoder(modified)
print(f'修改后解码: {decoded_modified}')
print(f'是否合法: {Chem.MolFromSmiles(decoded_modified) is not None}')
2.2 分子图表示
将分子视为图结构：原子为节点，化学键为边。这种表示天然适合图神经网络。
import torch
from torch_geometric.data import Data
def mol_to_graph(smiles):
"""将SMILES转换为PyTorch Geometric的图数据"""
mol = Chem.MolFromSmiles(smiles)
if mol is None:
return None
# 节点特征：原子属性
atom_features = []
for atom in mol.GetAtoms():
features = [
atom.GetAtomicNum(),
atom.GetDegree(),
atom.GetFormalCharge(),
int(atom.GetHybridization()),
int(atom.GetIsAromatic()),
atom.GetTotalNumHs(),
]
atom_features.append(features)
# 边索引和边特征
edge_index = []
edge_attr = []
for bond in mol.GetBonds():
i, j = bond.GetBeginAtomIdx(), bond.GetEndAtomIdx()
edge_index.extend([[i, j], [j, i]]) # 无向图
bond_features = [
int(bond.GetBondType()),
int(bond.GetIsConjugated()),
int(bond.IsInRing()),
]
edge_attr.extend([bond_features, bond_features])
x = torch.tensor(atom_features, dtype=torch.float)
edge_index = torch.tensor(edge_index, dtype=torch.long).t().contiguous()
edge_attr = torch.tensor(edge_attr, dtype=torch.float)
return Data(x=x, edge_index=edge_index, edge_attr=edge_attr)
# 测试
graph = mol_to_graph('CC(=O)Oc1ccccc1C(=O)O')
print(f'节点数: {graph.num_nodes}')
print(f'边数: {graph.num_edges}')
print(f'节点特征维度: {graph.x.shape}')
2.3 3D 构象表示
分子的三维构象决定了其与蛋白质靶点的结合方式。生成 3D 构象通常使用 RDKit 的距离几何算法或 ETKDG 方法。
def generate_3d_conformer(smiles):
"""生成分子的3D构象"""
mol = Chem.MolFromSmiles(smiles)
mol = Chem.AddHs(mol)
# 使用ETKDG算法生成3D构象
AllChem.EmbedMolecule(mol, AllChem.ETKDGv3())
AllChem.MMFFOptimizeMolecule(mol)
# 提取原子坐标
conf = mol.GetConformer()
positions = []
for i in range(mol.GetNumAtoms()):
pos = conf.GetAtomPosition(i)
positions.append([pos.x, pos.y, pos.z])
return np.array(positions), mol
coords, mol_3d = generate_3d_conformer('CC(=O)Oc1ccccc1C(=O)O')
print(f'原子坐标 shape: {coords.shape}')
print(f'第一个原子坐标: {coords[0]}')
3. 分子属性预测
3.1 消息传递神经网络（MPNN）
MPNN 是图神经网络的经典范式，通过在图上迭代传递和聚合邻居信息来学习节点（原子）表示，最终通过全局池化得到分子级别的表示。
import torch
import torch.nn as nn
from torch_geometric.nn import MessagePassing, global_mean_pool
from torch_geometric.data import Batch
class MPNNLayer(MessagePassing):
"""消息传递层"""
def __init__(self, node_dim, edge_dim, hidden_dim):
super().__init__(aggr='add')
self.message_fn = nn.Sequential(
nn.Linear(node_dim * 2 + edge_dim, hidden_dim),
nn.ReLU(),
nn.Linear(hidden_dim, hidden_dim)
)
self.update_fn = nn.GRU(hidden_dim, node_dim)
def forward(self, x, edge_index, edge_attr):
return self.propagate(edge_index, x=x, edge_attr=edge_attr)
def message(self, x_i, x_j, edge_attr):
return self.message_fn(torch.cat([x_i, x_j, edge_attr], dim=-1))
def update(self, aggr_out, x):
return self.update_fn(aggr_out.unsqueeze(0), x.unsqueeze(0)).squeeze(0)
class MolecularMPNN(nn.Module):
"""分子属性预测MPNN模型"""
def __init__(self, node_dim=6, edge_dim=3, hidden_dim=128, n_layers=3, n_tasks=1):
super().__init__()
self.node_encoder = nn.Linear(node_dim, hidden_dim)
self.layers = nn.ModuleList([
MPNNLayer(hidden_dim, edge_dim, hidden_dim) for _ in range(n_layers)
])
self.readout = nn.Sequential(
nn.Linear(hidden_dim, 64),
nn.ReLU(),
nn.Linear(64, n_tasks)
)
def forward(self, data):
x = self.node_encoder(data.x)
for layer in self.layers:
x = layer(x, data.edge_index, data.edge_attr) + x # 残差连接
x = global_mean_pool(x, data.batch)
return self.readout(x)
# 训练循环示例
model = MolecularMPNN(n_tasks=1)
optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)
# 假设已有train_loader包含分子图数据和标签
# for batch in train_loader:
# pred = model(batch).squeeze()
# loss = nn.MSELoss()(pred, batch.y)
# optimizer.zero_grad()
# loss.backward()
# optimizer.step()
3.2 图注意力网络（GAT）
GAT 通过注意力机制为不同邻居分配不同权重，能够自适应地捕捉原子间的相互作用重要性。
from torch_geometric.nn import GATv2Conv, global_mean_pool
class MolecularGAT(nn.Module):
def __init__(self, in_dim=6, hidden_dim=128, heads=4, n_layers=3, n_tasks=1):
super().__init__()
self.enc = nn.Linear(in_dim, hidden_dim)
self.convs = nn.ModuleList()
for _ in range(n_layers):
self.convs.append(GATv2Conv(hidden_dim, hidden_dim // heads, heads=heads))
self.head = nn.Linear(hidden_dim, n_tasks)
def forward(self, data):
x = torch.relu(self.enc(data.x))
for conv in self.convs:
x = torch.relu(conv(x, data.edge_index)) + x
x = global_mean_pool(x, data.batch)
return self.head(x)
4. 分子生成模型
4.1 分子 VAE（变分自编码器）
VAE 将分子编码到连续隐空间，然后从隐空间解码生成新分子。通过在隐空间中插值或优化，可以生成具有目标属性的分子。
class MolecularVAE(nn.Module):
"""基于SELFIES的分子VAE"""
def __init__(self, vocab_size, max_len, embed_dim=64, latent_dim=128):
super().__init__()
self.max_len = max_len
# 编码器
self.encoder = nn.GRU(embed_dim, 256, bidirectional=True, batch_first=True)
self.embed = nn.Embedding(vocab_size, embed_dim)
self.fc_mu = nn.Linear(512, latent_dim)
self.fc_logvar = nn.Linear(512, latent_dim)
# 解码器
self.latent_to_hidden = nn.Linear(latent_dim, 256)
self.decoder = nn.GRU(embed_dim, 256, batch_first=True)
self.fc_out = nn.Linear(256, vocab_size)
def encode(self, x):
embedded = self.embed(x)
_, h = self.encoder(embedded)
h = torch.cat([h[0], h[1]], dim=-1)
return self.fc_mu(h), self.fc_logvar(h)
def reparameterize(self, mu, logvar):
std = torch.exp(0.5 * logvar)
eps = torch.randn_like(std)
return mu + eps * std
def decode(self, z, max_len):
h = self.latent_to_hidden(z).unsqueeze(0)
input_token = torch.zeros(z.shape[0], 1, dtype=torch.long, device=z.device)
outputs = []
for _ in range(max_len):
embedded = self.embed(input_token)
out, h = self.decoder(embedded, h)
logits = self.fc_out(out)
outputs.append(logits)
input_token = logits.argmax(dim=-1)
return torch.cat(outputs, dim=1)
def forward(self, x):
mu, logvar = self.encode(x)
z = self.reparameterize(mu, logvar)
recon = self.decode(z, x.shape[1])
return recon, mu, logvar
def vae_loss(recon, x, mu, logvar, beta=0.5):
"""VAE损失 = 重建损失 + β * KL散度"""
recon_loss = nn.CrossEntropyLoss()(recon.view(-1, recon.shape[-1]), x.view(-1))
kl_loss = -0.5 * torch.sum(1 + logvar - mu.pow(2) - logvar.exp()) / x.shape[0]
return recon_loss + beta * kl_loss
4.2 分子扩散模型
扩散模型通过逐步去噪过程生成分子，在分子生成质量和多样性方面表现出色。
class MolecularDiffusion(nn.Module):
"""简化的分子属性条件扩散模型（在隐空间操作）"""
def __init__(self, latent_dim=128, time_dim=64, cond_dim=32):
super().__init()
self.time_embed = nn.Sequential(
nn.Linear(1, time_dim),
nn.SiLU(),
nn.Linear(time_dim, time_dim)
)
self.cond_embed = nn.Linear(cond_dim, time_dim)
self.net = nn.Sequential(
nn.Linear(latent_dim + time_dim * 2, 256),
nn.SiLU(),
nn.Linear(256, 256),
nn.SiLU(),
nn.Linear(256, latent_dim)
)
def forward(self, x, t, condition):
t_emb = self.time_embed(t.unsqueeze(-1).float())
c_emb = self.cond_embed(condition)
inp = torch.cat([x, t_emb, c_emb], dim=-1)
return self.net(inp)
def diffuse_sample(model, condition, latent_dim=128, n_steps=1000):
"""从扩散模型采样"""
device = next(model.parameters()).device
x = torch.randn(1, latent_dim, device=device)
cond = torch.tensor(condition, dtype=torch.float, device=device).unsqueeze(0)
for t in reversed(range(n_steps)):
t_tensor = torch.tensor([t / n_steps], device=device)
predicted_noise = model(x, t_tensor, cond)
# 简化的去噪步骤
alpha = 1 - t / n_steps
x = (x - predicted_noise * (1 - alpha)) / max(alpha, 1e-8)
if t > 0:
x += torch.randn_like(x) * 0.01
return x.detach()
4.3 基于强化学习的分子优化
将分子生成建模为马尔可夫决策过程，使用策略梯度方法优化生成的分子属性。
def compute_reward(smiles, target_props):
"""计算分子奖励分数"""
mol = Chem.MolFromSmiles(smiles)
if mol is None:
return -1.0 # 无效分子惩罚
reward = 0.0
# LogP奖励（药物样分子通常在-0.4到5.6之间）
logp = Descriptors.MolLogP(mol)
if -0.4 <= logp <= 5.6:
reward += 1.0
# 分子量惩罚
mw = Descriptors.MolWt(mol)
if 200 <= mw <= 500:
reward += 1.0
# 药物相似性（Lipinski规则）
if Descriptors.NumHDonors(mol) <= 5 and Descriptors.NumHAcceptors(mol) <= 10:
reward += 0.5
# 合成可及性（越低越容易合成）
# sa_score = calculate_sa_score(mol)
# reward += max(0, 10 - sa_score) / 10
return reward
5. 蛋白质结构预测
5.1 AlphaFold 核心思想
AlphaFold 通过 Evoformer（处理多序列比对 MSA 和配对表示）和 Structure Module（预测原子坐标）两个核心模块，实现了接近实验精度的蛋白质结构预测。
# 使用ColabFold API进行蛋白质结构预测（简化示例）
# 实际使用需要安装 colabfold 或访问 AlphaFold DB
def fetch_predicted_structure(protein_id):
"""从AlphaFold数据库获取预测结构"""
import urllib.request
url = f'https://alphafold.ebi.ac.uk/files/AF-{protein_id}-F1-model_v4.pdb'
pdb_path = f'{protein_id}.pdb'
urllib.request.urlretrieve(url, pdb_path)
return pdb_path
# 蛋白质序列编码（用于自定义模型）
def encode_protein_sequence(seq, max_len=1024):
"""将氨基酸序列编码为数值向量"""
aa_vocab = 'ACDEFGHIKLMNPQRSTVWY'
aa_to_idx = {aa: i+1 for i, aa in enumerate(aa_vocab)}
encoded = [aa_to_idx.get(aa, 0) for aa in seq[:max_len]]
# 填充到固定长度
encoded += [0] * (max_len - len(encoded))
return encoded
# 示例
seq = 'MKFLILLFNILCLFPVLAADNHGVSRELVDKGKLF'
encoded = encode_protein_sequence(seq)
print(f'序列长度: {len(seq)}, 编码后长度: {len(encoded)}')
5.2 ESMFold
ESMFold 是 Meta 开发的蛋白质结构预测模型，无需 MSA 即可预测结构，速度比 AlphaFold 快约 60 倍。
# ESMFold使用示例（需要安装esm库）
# pip install fair-esm
import esm
def predict_with_esmfold(sequence):
"""使用ESMFold预测蛋白质结构"""
model = esm.pretrained.esmfold_v1()
model = model.eval()
with torch.no_grad():
output = model.infer_pdb(sequence)
# 保存PDB文件
with open('predicted.pdb', 'w') as f:
f.write(output)
return 'predicted.pdb'
# 蛋白质语言模型嵌入（用于下游任务）
def get_protein_embedding(sequence):
"""使用ESM-2获取蛋白质序列嵌入"""
model, alphabet = esm.pretrained.esm2_t33_650M_UR50D()
batch_converter = alphabet.get_batch_converter()
data = [('protein', sequence)]
_, _, tokens = batch_converter(data)
with torch.no_grad():
results = model(tokens, repr_layers=[33])
# 取最后一层的表示，平均池化
embedding = results['representations'][33].mean(dim=1)
return embedding
# embed = get_protein_embedding('MKFLILLFNILCLFPVLAADNHGVSRELVD')
# print(f'蛋白质嵌入维度: {embed.shape}')
6. 药物-靶点相互作用预测
预测药物分子与蛋白质靶点的结合亲和力是药物发现的核心任务。
class DTIPredictor(nn.Module):
"""药物-靶点相互作用预测模型"""
def __init__(self, drug_dim=128, protein_dim=256, hidden_dim=128):
super().__init__()
# 药物编码器（MPNN）
self.drug_encoder = nn.Sequential(
nn.Linear(drug_dim, hidden_dim),
nn.ReLU(),
nn.Linear(hidden_dim, hidden_dim)
)
# 蛋白质编码器
self.protein_encoder = nn.Sequential(
nn.Linear(protein_dim, hidden_dim),
nn.ReLU(),
nn.Linear(hidden_dim, hidden_dim)
)
# 交互层
self.interaction = nn.Sequential(
nn.Linear(hidden_dim * 2, hidden_dim),
nn.ReLU(),
nn.Dropout(0.3),
nn.Linear(hidden_dim, 64),
nn.ReLU(),
nn.Linear(64, 1),
nn.Sigmoid()
)
def forward(self, drug_repr, protein_repr):
d = self.drug_encoder(drug_repr)
p = self.protein_encoder(protein_repr)
combined = torch.cat([d, p], dim=-1)
return self.interaction(combined)
# 分子指纹作为药物表示
def get_morgan_fingerprint(smiles, radius=2, n_bits=1024):
"""生成Morgan分子指纹"""
mol = Chem.MolFromSmiles(smiles)
fp = AllChem.GetMorganFingerprintAsBitVect(mol, radius, nBits=n_bits)
return np.array(fp)
# 使用示例
drug_fp = get_morgan_fingerprint('CC(=O)Oc1ccccc1C(=O)O')
print(f'分子指纹维度: {drug_fp.shape}')
print(f'非零位点数: {np.sum(drug_fp > 0)}')
7. 虚拟筛选与先导化合物优化
7.1 虚拟筛选流程
虚拟筛选从数百万化合物中快速缩小候选范围，分为基于结构的筛选（SBVS）和基于配体的筛选（LBVS）。
from rdkit.Chem import AllChem
from rdkit.DataStructs import TanimotoSimilarity
def virtual_screening(query_smiles, library_smiles, top_k=10):
"""基于配体相似性的虚拟筛选"""
query_mol = Chem.MolFromSmiles(query_smiles)
query_fp = AllChem.GetMorganFingerprintAsBitVect(query_mol, 2, nBits=2048)
results = []
for name, smi in library_smiles.items():
mol = Chem.MolFromSmiles(smi)
if mol is None:
continue
fp = AllChem.GetMorganFingerprintAsBitVect(mol, 2, nBits=2048)
similarity = TanimotoSimilarity(query_fp, fp)
results.append((name, smi, similarity))
results.sort(key=lambda x: x[2], reverse=True)
return results[:top_k]
# 示例：筛选与阿司匹林相似的化合物
library = {
'化合物A': 'CC(=O)Oc1ccccc1C(=O)O', # 阿司匹林本身
'化合物B': 'CC(=O)Nc1ccc(O)cc1', # 对乙酰氨基酚
'化合物C': 'CC(C)Cc1ccc(cc1)C(C)C(=O)O', # 布洛芬
'化合物D': 'c1ccc2c(c1)cc1ccccc1c2', # 芘
}
hits = virtual_screening('CC(=O)Oc1ccccc1C(=O)O', library)
for name, smi, sim in hits:
print(f'{name}: Tanimoto = {sim:.3f}')
7.2 先导化合物优化
通过 AI 生成分子变体，在保持核心骨架的同时优化 ADMET（吸收、分布、代谢、排泄、毒性）属性。
def lipinski_filter(smiles):
"""Lipinski五规则过滤"""
mol = Chem.MolFromSmiles(smiles)
if mol is None:
return False
return (Descriptors.MolWt(mol) <= 500 and
Descriptors.MolLogP(mol) <= 5 and
Descriptors.NumHDonors(mol) <= 5 and
Descriptors.NumHAcceptors(mol) <= 10)
def generate_analogs(base_smiles, n_samples=100):
"""基于SMILES扰动生成类似物"""
mol = Chem.MolFromSmiles(base_smiles)
analogs = set()
for _ in range(n_samples * 10):
# 随机原子替换和键修改
new_mol = Chem.RWMol(mol)
# 简化：对SMILES字符串做微小修改
chars = list(base_smiles)
idx = np.random.randint(len(chars))
replacement = np.random.choice(['C', 'N', 'O', 'F', 'Cl'])
chars[idx] = replacement
new_smiles = ''.join(chars)
new_mol_check = Chem.MolFromSmiles(new_smiles)
if new_mol_check is not None and lipinski_filter(new_smiles):
analogs.add(new_smiles)
if len(analogs) >= n_samples:
break
return list(analogs)
# analogs = generate_analogs('CC(=O)Oc1ccccc1C(=O)O', n_samples=10)
# print(f'生成了 {len(analogs)} 个合法类似物')
8. 分子动力学模拟
分子动力学（MD）模拟追踪分子体系随时间的运动轨迹，用于评估蛋白质-配体复合物的稳定性。AI 正在加速 MD 模拟中的力场计算和采样。
# MD模拟概念性代码（使用OpenMM的简化接口）
# pip install openmm
def run_md_simulation_concept():
"""MD模拟的核心流程（概念性展示）"""
# 1. 体系准备
steps = {
'加载拓扑': '从PDB文件读取原子坐标和拓扑',
'添加溶剂': '将分子置于水盒子中',
'添加离子': '平衡体系电荷',
'能量最小化': '消除原子间的空间冲突',
'NVT平衡': '恒温平衡，使体系达到目标温度',
'NPT平衡': '恒温恒压平衡，调整体系密度',
'生产运行': '收集轨迹数据用于分析',
}
# 2. AI加速的力场（概念）
# 传统力场：基于预定义的数学函数
# ML力场：从量子力学数据学习势能面
# 代表模型：ANI, SchNet, MACE
# 3. 关键分析指标
metrics = {
'RMSD': '均方根偏差 - 评估结构变化',
'RMSF': '均方根波动 - 评估残基柔性',
'Rg': '回旋半径 - 评估蛋白质紧密度',
'氢键数': '评估结合稳定性',
'MM/PBSA': '结合自由能计算',
}
return steps, metrics
steps, metrics = run_md_simulation_concept()
print('MD模拟流程:')
for k, v in steps.items():
print(f' {k}: {v}')
AI 加速 MD 的关键方向：
- ML 力场：用图神经网络替代传统力场，精度接近量子力学，速度接近经典力场
- 增强采样：用 AI 学习集体变量，加速稀有事件的采样
- 粗粒化模型：通过深度学习将全原子模型映射为粗粒化表示
9. 临床试验 AI 辅助
AI 在临床试验阶段的应用日益重要：
def clinical_trial_ai_applications():
"""AI在临床试验中的应用场景"""
applications = {
'患者招募': {
'技术': 'NLP分析电子病历，匹配入排标准',
'效果': '招募速度提升30-50%',
},
'剂量优化': {
'技术': '贝叶斯自适应设计，PK/PD建模',
'效果': '减少所需受试者数量，缩短试验周期',
},
'终点预测': {
'技术': '生存分析模型，预测治疗反应',
'效果': '更准确的样本量估计',
},
'安全性监测': {
'技术': '实时不良事件信号检测',
'效果': '早期发现安全性问题',
},
'数字孪生': {
'技术': '基于历史数据生成虚拟对照组',
'效果': '减少安慰剂组人数',
},
}
return applications
# 患者分层示例
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import cross_val_score
def patient_stratification(features, labels):
"""基于生物标志物的患者分层"""
clf = RandomForestClassifier(n_estimators=100, random_state=42)
scores = cross_val_score(clf, features, labels, cv=5, scoring='roc_auc')
print(f'患者分层 AUC: {scores.mean():.3f} ± {scores.std():.3f}')
return clf
10. 实战案例：分子性质预测与生成
以下构建一个完整的分子性质预测与生成流水线。
import torch
import torch.nn as nn
import numpy as np
from rdkit import Chem
from rdkit.Chem import AllChem, Descriptors
from torch_geometric.data import Data, DataLoader
from torch_geometric.nn import GCNConv, global_mean_pool
# === 第一步：数据准备 ===
def prepare_molecular_dataset(smiles_list, property_fn):
"""准备分子图数据集"""
dataset = []
for smi in smiles_list:
mol = Chem.MolFromSmiles(smi)
if mol is None:
continue
graph = mol_to_graph(smi) # 使用前面定义的函数
if graph is None:
continue
graph.y = torch.tensor([property_fn(mol)], dtype=torch.float)
dataset.append(graph)
return dataset
# 生成示例数据
np.random.seed(42)
sample_smiles = [
'CC(=O)Oc1ccccc1C(=O)O', 'CC(=O)Nc1ccc(O)cc1',
'CC(C)Cc1ccc(cc1)C(C)C(=O)O', 'Cn1c(=O)c2c(ncn2C)n(C)c1=O',
'c1ccc(-c2ccccc2)cc1', 'CC(=O)c1ccc(cc1)C(=O)C',
'O=c1ccc2ccccc2o1', 'c1ccc2c(c1)[nH]c1ccccc12',
'CC12CCC3C(C1CCC2O)CCC4=CC(=O)CCC34C', 'CC(=O)OC1CC2CCC3C(CCC4(C)C3CCC4(O)C(=O)CO)C2(C)CC1O',
]
dataset = prepare_molecular_dataset(sample_smiles, Descriptors.MolLogP)
print(f'数据集大小: {len(dataset)}')
print(f'样本: {dataset[0]}')
# === 第二步：模型定义 ===
class MolecularGNN(nn.Module):
def __init__(self, in_dim, hidden_dim=64, out_dim=1):
super().__init__()
self.conv1 = GCNConv(in_dim, hidden_dim)
self.conv2 = GCNConv(hidden_dim, hidden_dim)
self.conv3 = GCNConv(hidden_dim, hidden_dim)
self.fc = nn.Linear(hidden_dim, out_dim)
self.dropout = nn.Dropout(0.2)
def forward(self, data):
x, edge_index, batch = data.x, data.edge_index, data.batch
x = torch.relu(self.conv1(x, edge_index))
x = self.dropout(x)
x = torch.relu(self.conv2(x, edge_index))
x = self.dropout(x)
x = torch.relu(self.conv3(x, edge_index))
x = global_mean_pool(x, batch)
return self.fc(x)
# === 第三步：训练 ===
def train_model(dataset, epochs=100):
train_size = int(0.8 * len(dataset))
train_data = dataset[:train_size]
test_data = dataset[train_size:]
train_loader = DataLoader(train_data, batch_size=16, shuffle=True)
test_loader = DataLoader(test_data, batch_size=16)
model = MolecularGNN(in_dim=dataset[0].x.shape[1])
optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)
criterion = nn.MSELoss()
for epoch in range(epochs):
model.train()
total_loss = 0
for batch in train_loader:
pred = model(batch).squeeze()
loss = criterion(pred, batch.y)
optimizer.zero_grad()
loss.backward()
optimizer.step()
total_loss += loss.item()
if (epoch + 1) % 20 == 0:
model.eval()
test_loss = 0
with torch.no_grad():
for batch in test_loader:
pred = model(batch).squeeze()
test_loss += criterion(pred, batch.y).item()
print(f'Epoch {epoch+1}, Train Loss: {total_loss/len(train_loader):.4f}, '
f'Test Loss: {test_loss/len(test_loader):.4f}')
return model
# model = train_model(dataset, epochs=100)
# === 第四步：分子生成 ===
def generate_novel_molecules(base_smiles, model, n_generate=10):
"""基于梯度引导的分子生成"""
generated = []
base_mol = Chem.MolFromSmiles(base_smiles)
# 简化：通过随机SMILES变异生成候选
chars = list('CNOSF()=c123456789')
attempts = 0
while len(generated) < n_generate and attempts < 1000:
attempts += 1
smiles_chars = list(Chem.MolToSmiles(base_mol))
# 随机修改1-2个字符
n_changes = np.random.randint(1, 3)
for _ in range(n_changes):
idx = np.random.randint(len(smiles_chars))
smiles_chars[idx] = np.random.choice(chars)
new_smiles = ''.join(smiles_chars)
new_mol = Chem.MolFromSmiles(new_smiles)
if new_mol is not None:
generated.append((new_smiles, Descriptors.MolLogP(new_mol)))
return generated
# novel = generate_novel_molecules('CC(=O)Oc1ccccc1C(=O)O', model)
# for smi, logp in novel[:5]:
# print(f'{smi} -> LogP: {logp:.2f}')
11. 挑战与前沿方向
当前挑战
- 数据稀缺性：高质量标注的药物活性数据极为有限，尤其针对新靶点
- 化学空间巨大：类药物分子的化学空间约有 10^60 种可能，穷举搜索不可行
- 多目标优化：需要同时优化活性、选择性、ADMET 属性、合成可行性等多个目标
- 实验验证鸿沟：计算预测与湿实验结果之间仍存在显著差距
- 可解释性：黑盒模型的预测缺乏化学直觉支撑，难以获得药物化学家信任
前沿方向
- 基础模型：如 MolFormer、Uni-Mol 等在大规模分子数据上预训练的通用模型
- 多模态学习：融合分子图、SMILES、3D 构象、蛋白质序列等多种模态
- 因果推理：从相关性建模转向因果机制理解，提升模型泛化能力
- 主动学习：智能选择最有信息量的实验进行验证，最大化实验效率
- 量子机器学习：利用量子计算加速分子模拟和优化
- 闭环实验：AI 设计 → 机器人合成 → 自动检测 → 反馈优化的全自动流程
AI 药物发现正处于从概念验证到产业化落地的关键转折期。掌握分子表示学习、属性预测和生成模型三大核心技术，理解从计算到实验的完整闭环，是进入这一领域的基础。随着基础模型和自动化实验平台的发展，AI 辅助药物发现的效率还将持续提升。
