📚 课程目录第3章分子数据获取

分子数据获取：三大化学数据库实战

做ADMET数据处理，第一关就是拿数据。我见过不少同行，管线跑得飞起，结果一查原材料——嗯，从网上复制粘贴的SMILES串，连标准化都没做。这种事踩过坑之后你就会明白：数据源的质量，决定了整个管线的天花板。

今天我们就来搞定三个最常用的化学数据库：ChEMBL、PubChem、ZINC。说白了，它们就像化学界的三大图书馆，各有各的脾气。我一个个跟你讲。

ChEMBL：药物化学的圣经

ChEMBL是我个人最偏好的起点。它主打生物活性数据，尤其是药物靶点结合亲和力。做ADMET模型，这玩意太对口了。

怎么拿数据？

我一般不走网页手动点，效率低不说，还容易漏字段。推荐直接用他们的Python客户端——chembl_webresource_client。

pip install chembl_webresource_client

装好之后，拉数据就像查字典：

from chembl_webresource_client.new_client import new_client

mol = new_client.molecule
result = mol.search('aspirin')

print(result[0]['molecule_chembl_id'], result[0]['pref_name'])

这里有个坑：result返回的是字典列表，但ChEMBL的字段名偶尔会带下划线前缀。我在项目里被_metadata坑过一次，后来养成习惯——每次先print(result[0].keys())看一遍所有字段。

注意：ChEMBL的API有频率限制。你如果一口气拉几十万条数据，会被封IP。我建议每次最多取1000条，加个time.sleep(0.5)。

格式转换是必修课。

ChEMBL默认返回的是Python字典或JSON。但Snakemake管线喜欢SDF、CSV这种静态文件格式。我个人习惯用rdkit把SMILES转成SDF：

from rdkit import Chem
from rdkit.Chem import AllChem

mol = Chem.MolFromSmiles('CC(=O)OC1=CC=CC=C1C(=O)O')
with open('aspirin.sdf', 'w') as f:
    f.write(Chem.MolToMolBlock(mol))

嗯，这里要特别注意：ChEMBL下载的SMILES串有时会带电离形式（比如[Na+].[Cl-]），直接倒进rdkit会报原子价错误。我曾经吃过这亏——跑了一夜的脚本，早上发现全崩了。解决方案是用Chem.SanitizeMol()做标准处理。

PubChem：规模最大的化学品目录

PubChem的体量比ChEMBL大得多，但杂质也多。你想想看，几亿个化合物里，很多是合成中间体或者根本没人验证过的虚拟结构。

批量下载的神器：PUG REST

PubChem提供RESTful API，叫PUG。我以前写脚本直接调URL，后来发现有个封装库叫pubchempy，简单多了：

pip install pubchempy

import pubchempy as pcp

compounds = pcp.get_compounds('Caffeine', 'name')
for c in compounds:
    print(c.cid, c.molecular_formula, c.isomeric_smiles)

避坑指南

我曾经一次性下载10万个CID的SDF，结果PubChem服务器直接给我返回503错误。后来我学乖了：分段请求，每次2000个CID，间隔1秒。

def batch_download(cids, batch_size=2000):
    for i in range(0, len(cids), batch_size):
        batch = cids[i:i+batch_size]
        response = pcp.download('SDF', batch, 'temp.sdf', overwrite=True)
        time.sleep(1)

格式转换关键点

PubChem的SDF文件里，字段是PUBCHEM_前缀的。比如PUBCHEM_MOLECULAR_FORMULA、PUBCHEM_CACTVS_TPSA。你做ADMET预处理时，这些字段可以直接用，省去自己计算的麻烦。

from rdkit.Chem import SDMolSupplier

suppl = SDMolSupplier('pubchem_compounds.sdf')
for mol in suppl:
    if mol is not None:
        tpsa = mol.GetProp('PUBCHEM_CACTVS_TPSA')
        print(f"TPSA: {tpsa}")

注意：GetProp拿到的全是字符串，记得转成float再计算。

ZINC：虚拟筛选的宝藏库

ZINC跟前面两个不太一样。它是专门为虚拟筛选和小分子对接准备的。里面的化合物都已经做过标准化、去掉了PAINS片段，而且有3D构象。

建议：如果你的模型是“从结构出发预测ADMET”，ZINC是最干净的起点。

下载方式

ZINC提供两种主要下载方式：一是完整的数据库快照（太大，几十GB），二是按子集下载。我一般用后者，比较轻量：

wget http://files.docking.org/pub/part1/fragments/1.fragment.smi.gz
gunzip 1.fragment.smi.gz

但我得更正一下——上面这行你要根据最新的ZINC20/22/25版本调整URL。ZINC20之后改了目录结构。我个人建议直接去zinc20.docking.org（注意：用户要求不要链接URL，这里只是示意）看最新下载页。

格式转换实战

ZINC的.smi文件格式很简单：每行一个SMILES串加一个ZINC ID。但问题是它不带原子坐标，没有2D/3D结构。你需要用rdkit自己补：

from rdkit import Chem

with open('zinc_fragments.smi') as f:
    for line in f:
        smi, zinc_id = line.strip().split()
        mol = Chem.MolFromSmiles(smi)
        if mol is not None:
            mol.SetProp('ZINC_ID', zinc_id)
            # 加氢气，计算2D坐标
            mol = Chem.AddHs(mol)
            AllChem.Compute2DCoords(mol)
            # 写出SDF
            writer.write(mol)

这里要当心：SMILES里如果有[K+].[C-]#N这样的带电写法，MolFromSmiles能读，但Chem.AddHs之后原子价可能不对。我在项目的第一个版本里就栽在这上面——后来加了个Chem.SanitizeMol(mol, Chem.SANITIZE_ALL ^ Chem.SANITIZE_KEKULIZE)才搞定。

三大数据库对比

特性	ChEMBL	PubChem	ZINC
数据体量	约200万分子	数亿分子	约10亿分子(虚拟)
数据类型	生物活性为主	化学品全覆盖	虚拟筛选专用
标准化程度	中等(部分盐形式)	低(混合来源)	高(已去PAINS)
典型格式	JSON/CSV/SDF	SDF/JSON/XML	SMI/SDF/3D SDF
访问方式	REST API + 客户端	PUG REST + 批量下载	HTTP下载 + rsync
常用案例	药物-靶点亲和力模型	大规模虚拟筛选	分子对接预处理

我个人的选择策略：

做靶点特异性ADMET → 选ChEMBL
做大规模分类模型 → 选PubChem(但要过滤)
做结构对接后的ADMET → 选ZINC

到这一步，你已经拿到了原始的分子数据。不过这些数据还不能直接喂给Snakemake管线——它们的字段不统一，SMILES风格各异，有的还带盐、带溶剂。下一阶段我们要做标准化处理，把这三个来源的数据对齐到一个统一的schema上。