rdkit作为一个强大的化学信息学工具包，不仅能够处理分子结构数据，还提供了丰富的分子可视化功能。在药物发现和化学性质研究中，识别和可视化分子的极性区域至关重要，因为它直接影响分子的溶解度、渗透性以及与生物靶点的相互作用。本教程将详细介绍几种在rdkit中实现分子极性区域可视化的方法，从简单的原子高亮到更复杂的相似性图，帮助用户根据需求选择最合适的策略。
一、基于Gasteiger电荷的原子高亮显示
最初，用户可能倾向于通过计算原子上的部分电荷来识别极性原子。Gasteiger电荷模型是一种常用的半经验方法，用于估计分子中原子的部分电荷。通过识别带有显著负电荷（或正电荷）的原子，可以尝试将其高亮显示。
1.1 实现方法
以下代码演示了如何计算分子的Gasteiger电荷，并高亮显示带有负电荷的原子：
from rdkit import Chem
from rdkit.Chem import Draw
from rdkit.Chem import AllChem
def highlight_psa_atoms_by_gasteiger(mol):
"""
根据Gasteiger电荷高亮显示具有负电荷的原子。
"""
AllChem.ComputeGasteigerCharges(mol)
# 筛选出Gasteiger电荷小于0的原子
psa_atoms = [atom.GetIdx() for atom in mol.GetAtoms() if atom.GetDoubleProp("_GasteigerCharge") < 0]
# 设置高亮颜色（红色）
highlight_style = {atom_id: (1, 0, 0) for atom_id in psa_atoms}
return highlight_style
# 示例分子：阿司匹林
smiles = "CC(=O)OC1=CC=CC=C1C(O)=O"
mol = Chem.MolFromSmiles(smiles)
# 获取高亮样式
highlight_style = highlight_psa_atoms_by_gasteiger(mol)
# 绘制分子并高亮显示
img = Draw.MolToImage(mol, size=(300, 300), highlightAtoms=highlight_style, wedgeBonds=True, kekulize=True, wedgeLineWidth=2)
img1.2 注意事项与局限性
这种方法虽然直观，但存在一定的局限性。Gasteiger电荷计算可能会在芳香环等区域产生小的部分电荷，导致这些区域也被错误地高亮显示，而它们通常不被认为是分子极性表面积（Polar Surface Area, PSA）的贡献者。PSA通常特指由氧和氮原子（以及硫和磷，取决于定义）承载的极性氢原子所构成的表面积。因此，简单地基于Gasteiger电荷进行高亮，可能无法准确反映化学上通常理解的极性区域。
二、基于总极性表面积（TPSA）贡献度的可视化
为了更准确地识别分子的极性区域，RDKit提供了计算总极性表面积（TPSA）的功能。TPSA是一个重要的描述符，它量化了分子中极性原子（通常是氧和氮）的表面积。RDKit能够计算每个原子对TPSA的贡献，这为精确可视化极性区域提供了基础。
2.1 实现方法
RDKit的rdMolDescriptors._CalcTPSAContribs函数可以计算每个原子对TPSA的贡献值。贡献值大于0的原子即为TPSA的贡献者，通常是极性原子。
from rdkit import Chem from rdkit.Chem import Draw from rdkit.Chem.Draw import rdMolDraw2D from rdkit.Chem import rdMolDescriptors from IPython.display import Image # 示例分子：阿司匹林 smiles = "CC(=O)OC1=CC=CC=C1C(O)=O" mol = Chem.MolFromSmiles(smiles) # 计算每个原子对TPSA的贡献 # includeSandP=True 参数可以控制是否将硫和磷原子也纳入TPSA计算 tpsa_contribs = rdMolDescriptors._CalcTPSAContribs(mol, includeSandP=True) # 筛选出对TPSA有贡献的原子索引 highlight_atoms = [i for i, contrib in enumerate(tpsa_contribs) if contrib > 0] # 创建一个绘图对象，用于生成PNG图像 drawer = rdMolDraw2D.MolDraw2DCairo(300, 300) # 绘制分子并高亮显示贡献原子 drawer.DrawMolecule(mol, highlightAtoms=highlight_atoms) drawer.FinishDrawing() # 获取PNG数据并显示 png_data = drawer.GetDrawingText() Image(png_data)
2.2 参数说明
- rdMolDescriptors._CalcTPSAContribs(mol, includeSandP=True): 这个函数返回一个列表，其中每个元素对应分子中一个原子对TPSA的贡献值。includeSandP参数用于指定是否将硫（S）和磷（P）原子也视为极性原子并纳入TPSA计算。默认情况下，RDKit可能不包含它们。
这种方法能够更准确地识别出对极性表面积有实际贡献的原子，避免了Gasteiger电荷方法中可能出现的误判。
三、使用相似性图（Similarity Maps）进行渐变色可视化
除了简单的原子高亮，RDKit还提供了更高级的可视化工具——相似性图（Similarity Maps），可以将原子权重映射到分子结构上，通过渐变色来表示不同区域的极性强度，从而提供更直观、更细致的极性分布视图。
3.1 实现方法
SimilarityMaps.GetSimilarityMapFromWeights函数允许用户将一个权重列表（例如TPSA贡献值）应用于分子，并生成一个带有颜色渐变的图像。
import numpy as np
from rdkit import Chem
from rdkit.Chem import AllChem, Draw, rdMolDescriptors
from rdkit.Chem.Draw import SimilarityMaps
from matplotlib import pyplot as plt # 导入pyplot用于显示图像
# 示例分子
smiles = "CCNC(=O)NC1=NC2=CC=C(C=C2S1)C(=O)NCCS"
mol = Chem.MolFromSmiles(smiles)
# 计算每个原子对TPSA的贡献
# includeSandP=True 参数可以控制是否将硫和磷原子也纳入TPSA计算
tpsa_contribs = rdMolDescriptors._CalcTPSAContribs(mol, includeSandP=True)
# 生成相似性图
# weights: 每个原子的权重，这里使用TPSA贡献值
# colorMap: 颜色映射方案，例如 'bwr' (蓝白红) 是一种常见的发散型色图
# contourLines: 等高线数量，用于在图上显示不同权重的分界线
fig = SimilarityMaps.GetSimilarityMapFromWeights(
mol,
size=(400, 400),
weights=tpsa_contribs,
colorMap='bwr', # 可以尝试其他色图，如 'viridis', 'plasma', 'jet' 等
contourLines=10 # 调整等高线数量以获得最佳视觉效果
)
# 保存图像或直接显示
fig.savefig('tpsa_similarity_map.png', bbox_inches='tight')
plt.show() # 在Jupyter Notebook或脚本中显示图像3.2 参数说明与最佳实践
- weights: 这是一个与分子中原子数量相同的数值列表，每个值代表对应原子的权重。在本例中，我们使用TPSA贡献值。
- colorMap: 颜色映射方案。对于表示极性强度等连续数据，推荐使用发散型色图（如bwr，蓝代表低值，红代表高值，白色代表中间值）或顺序型色图（如viridis, plasma等）。matplotlib库提供了丰富的色图选择，可以在其官方文档中查看。
- contourLines: 等高线的数量。等高线可以帮助用户更清晰地识别不同极性强度区域的边界。
- 选择合适的色图：对于极性分布，通常希望能够区分“非极性”到“极性”的连续变化。bwr（蓝白红）是一个很好的选择，它可以清晰地表示正负或低高两端。
四、总结与建议
本教程介绍了RDKit中可视化分子极性区域的三种主要方法：
- 基于Gasteiger电荷的原子高亮：简单易行，但可能不够精确，容易将非极性区域误判为极性。适用于快速初步的定性分析。
- 基于TPSA贡献度的原子高亮：更准确地识别对总极性表面积有贡献的原子，是推荐的精确识别极性原子的方法。
- 使用相似性图进行渐变色可视化：提供最直观、最细致的极性强度分布视图，通过颜色渐变而非简单的点或线来表示极性。适用于需要深入分析极性分布模式的场景。
在实际应用中，推荐优先使用基于TPSA贡献度的方法来精确识别极性原子，或使用相似性图来获得更丰富的极性分布信息。虽然RDKit在可视化美观度上可能不如一些商业软件（如OpenEye）那样精致，但其功能强大且高度可定制，足以满足大多数科学研究和药物设计中的需求。通过灵活运用这些工具，用户可以有效地在RDKit中探索和展示分子的极性特征。
