化学信息工程入门：通过PubChemPy搜索和RDKit分析胆汁酸化学信息

前言

文本由Putuma P. Gqamana博士撰写的简短自学教程。主要通过使用PubChem开源数据库，RDKit库等探索使用Python提高工作效率。

什么是PubChem和PubChemPy？

PubChemPy是一个Python模块，用于与PubChem进行交互。PubChem是一个由美国国家生物技术信息中心（NCBI）维护的大规模化学信息库。该库通过Web接口提供按名称、子结构和相似性搜索化学信息的功能。PubChem是全球最大的此类集合，并且是连接其他专门化学数据库的交叉节点，例如人类代谢物数据库（HMDB）、NIST化学手册等。每个化学标识符都可以检索到化学和物理属性、光谱库、生物活性、安全性和毒性信息、专利、文献引用等内容，支持机器学习所需的格式。

PubChemPy的主要功能包括：

• 化学标准化： 它可以标准化同一化学实体的不同表示形式。

• 化学文件格式转换： 支持多种化学文件格式，并能够在它们之间进行转换。

• 化学属性的提取和绘制： 它可以生成化学物质的可视化图像，并检索各种化学属性。

• 与pandas DataFrame的集成： 它可以使用pandas DataFrame构建属性表。

• 支持二维和三维坐标： 它可以为化合物生成二维和三维坐标。

• 检索IUPAC名称、商品名和同义词： 它可以检索给定化合物的IUPAC系统名称、商品名以及所有已知的同义词。

PubChemPy完全依赖于PubChem数据库和通过其PUG REST Web服务提供的化学工具包。这意味着它需要持续的互联网连接，这可能导致某些任务的执行速度较慢，尤其是在本地计算机上执行时。然而，它也意味着您可以利用PubChem数据库和化学工具包的广泛资源，通过其他工具包如RDKit进一步探究化学问题，正如本文中所展示的那样。

什么是RDKit？

RDKit是一个开源的化学信息学工具包，采用C++和Python编写，提供了一系列化学信息学和机器学习软件功能。

RDKit的一些关键特点：

• C++中的核心数据结构和算法：RDKit提供了一套强大的数据结构，用于表示分子、反应及相关化学概念。此外，它还包括各种算法，用于操控和分析这些数据结构。

• 通过Boost.Python生成的Python 3.x包装器：RDKit为其C++库提供了Python绑定，使得用户能够使用Python与RDKit进行交互。

• 通过SWIG生成的Java和C#包装器：除了Python，RDKit还为Java和C#提供了绑定，扩展了其在不同编程语言中的可用性。

• 最重要功能的JavaScript包装器：RDKit为一些最重要的功能提供了JavaScript包装器，使其可以在基于Web的应用程序中使用。

• 2D和3D分子操作：RDKit支持多种2D和3D分子操作，包括生成2D和3D坐标、分子对齐以及构象生成。

• 用于机器学习的描述符生成：RDKit能够生成多种分子描述符，这些描述符通常作为化学信息学中的机器学习模型的输入。

• PostgreSQL的分子数据库卡片：RDKit提供了一种卡片，可以让PostgreSQL理解化学结构并直接在数据库中执行化学信息学操作。

• 适用于商业的BSD许可证：RDKit在商业环境中友好使用，支持Mac、Windows和Linux操作系统。它在学术界和工业界广泛用于化学信息学领域的研究和开发。

例子：结合PubChemPy和RDKit的应用

在这个教程中，使用PubChemPy和RDKit来分析胆汁酸，这些物质在生物化学和临床科学中非常重要。以下是如何进行操作的过程：

• PubChemPy：这个Python模块用于通过化学名称等标识符（如胆汁酸的名称）在PubChem中搜索化学数据。搜索结果被存储在CSV文件中，包含如下信息：分子名称同义词SMILES（简化分子输入线性结构式）InChI（国际化学标识符）键分子量（MW）拓扑表面积（TPSA）亲脂性（XLogP）分子式

• 数据汇总：将搜索到的数据整理并汇总到一个本地的Excel电子表格中。这个新表格包含了更全面的分子信息，可用于进一步分析。

• 2D可视化：使用RDKit生成2D分子结构，并将其渲染在pandas的DataFrame中。

• 结构分析：RDKit可以用来分析结构基序（特定的原子或功能基团的排列）和计算不同胆汁酸分子间的Tanimoto相似性（用于衡量分子结构相似性的一种方法）。

胆汁酸在生物化学中非常重要，因为它们参与脂肪和脂溶性维生素的消化吸收，并在胆固醇代谢中起着调节作用。胆汁酸对肝脏功能的健康也至关重要。由于它们在临床上的重要性，特别是在与肝脏、代谢和胆固醇相关的疾病中，研究其化学性质和结构是非常有意义的。通过RDKit，可以识别胆汁酸的结构特征，并使用相似性度量方法（如Tanimoto相似性）探索它们与其他化学物质的关系。

什么是胆汁酸？

胆汁酸是主要存在于哺乳动物胆汁中的类固醇酸。它们在肝脏中合成，传统上以在消化和吸收膳食脂肪中的作用而闻名。由于其多种生理功能，胆汁酸也作为多种疾病的生物标志物引起了广泛关注，这些功能包括调节葡萄糖和脂质代谢、影响肠道微生物群组成以及作为信号分子发挥作用。

胆汁酸是胆固醇代谢的下游产物。因此，初级胆汁酸在早期阶段的代表性结构是类似胆固醇的分子，末端为一个羧酸尾部（如下所示）。次级胆汁酸由于小肠内肠道微生物群的活动，在结构上有更多的修饰，有时甚至在两个端点都有修饰。

最近的研究强调了胆汁酸在肠道微生物群相互作用中的作用，这对炎症性肠病（IBD）具有重要意义。胆汁酸组成的变化可以影响肠道炎症和屏障功能，使其成为IBD诊断和治疗反应的潜在生物标志物。此外，某些肠道微生物群产生的新型和次级胆汁酸，如最近对日本百岁老人肠道微生物群的研究所建议，可能具有抗衰老的治疗特性。

总之，胆汁酸作为多种疾病的生物标志物具有潜力，能够为疾病机制和潜在的治疗靶点提供深刻的见解。它们在代谢、肝脏功能和肠道健康中的多样性作用凸显了其在研究和临床诊断中的重要价值。

本次教程主要内容

• 学习和测试PubChemPy搜索功能，通过检索胆汁酸列表中的RefMet名称来获取化学信息数据表。以前，我使用pubchempy.get_compounds()函数逐个调用，逐个搜索标识符，然后通过列表推导式为每个特性或属性列组装记录的数据框。这一次，我利用了该函数的内置能力，通过调用关键字参数“as_dataframe = True”来直接返回数据框（是的，我阅读了API文档）。之后，我按化合物循环，将数据框添加到空列表中，最终将各个化合物的数据框合并成一个大的数据框。我还在完整的数据框中添加了分子（“mol”）对象，以便使用RDKit进行进一步的化学信息学分析。

• 应用学习（第1点），对以“Tauro…”开头的名称的子集进行操作。在上述操作基础上，我有条件地选择了列表中所有的牛磺胆酸类（即14个胆汁酸），这些是次级胆汁酸，通常会发生重度修饰，具体内容如下。

• 进行描述性统计分析，检查这些分子的数值分布。我使用了pandas函数<describe()>对两个数据框（子集和总集）中的所有数值数据进行描述性统计。此外，我还使用了matplotlib.pyplot绘制了直方图，以查看两个集合中的数值属性分布。例如，高度集中的数据集通常分布较紧，而全面的数据集则会显示出更广泛的分布。

• 提取SMARTS字符串，以定义并表示两个数据集中的平均结构中的相似性模式。SMARTS字符串是一种用于表示一组相对相似分子核心结构模式的有用方式。在本例中，SMARTS可能类似于结构构建块——胆烷结构，之后在胆固醇/胆汁酸代谢和排泄过程中进一步修饰。

• 使用从大数据集中通过目视检查选择的参考分子（mol对象）进行Tanimoto相似度分析。Tanimoto相似度是化学信息学中用来评估分子之间相似性的度量，类似于Jaccard相似性。它通过参考分子和目标分子之间的分子碎片/成分的重叠来评估相似性。相似度值范围从0.0（无相似性）到1.0（完全相同）。

• 在大数据框中添加一个包含相似度值的额外列。

• 最后，将数据框导出为Excel电子表格，以便必要时进行进一步分析。

我在本教程中使用的代码是用Python编写的，并在Jupyter Notebooks中呈现。我通常会为将机器学习输出保存为Excel电子表格留出空间，这已成为一种习惯。这样做可能是一个有用的表现方式，并且对于不习惯Jupyter Notebooks输出的读者或合作者来说更为方便。请参阅下面的代码块：

# 导入依赖库

import os
import requests
from io import BytesIO
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import pandas as pd
from pandas import *
import pubchempy as pubchem
from pubchempy import *
import rdkit as rdkit
from rdkit import Chem
from rdkit.Chem import PandasTools
from rdkit.Chem import rdFMCS
from rdkit.Chem import AllChem
from rdkit.Chem import Draw
from rdkit.Chem.Descriptors import ExactMolWt 
from rdkit.Chem.Draw import IPythonConsole
IPythonConsole.ipython_useSVG = True
from IPython.display import IFrame
from PIL import Image
import mols2grid
import py3Dmol

# 上传数据（我将csv文件上传到相同的文件夹中）

df = read_csv("bile_acids_nomenclature.csv")
acids = df.dropna().copy()

# 从PubChem调用中选择列
"""
可以选择的列有很多。我提前知道从单一标识符搜索的实践中，哪些属性是需要选择的。
"""

use_calls = ['canonical_smiles','h_bond_acceptor_count','h_bond_donor_count','inchi','inchikey','iupac_name','molecular_formula','molecular_weight',
       'tpsa','xlogp']

# 只选择Tauros
"""
我已经在文本文件中看到，大多数牛磺酸结合酸以"Tauro"开头。
正则表达式可能是更好的搜索方法，具体取决于可用的命名规则。
在这个例子中，输出预计是一个简单的1D列表，
包含以字符串形式表示的RefMet名称，因为后续将使用的get_compounds()函数
仅接受字符串作为标识符。
"""

tauros = acids[acids["RefMet Name"].str.startswith("Tauro")]["RefMet Name"]

# 准备一个新数据框，仅包含所选的标识符记录，例如，牛磺胆酸

"""
策略是设置一个空列表，然后通过名称循环pubchem.get_compounds()搜索。
之后，我将数据框附加到空列表中。最后，我将各个分子数据框合并成一个大数据框。

警告：需要注意的5个问题：
1. 单次搜索可能返回“CID”作为索引列，这会修改后续循环并导致其崩溃，
因为与选择列表不匹配。
2. 但是，可能会发生相反的情况，即不会返回CID。会造成同样的不便。
3. 如果选择列表中包括“cid”，则会出现Key Errors（“cid”）和不一致的索引问题。
4. 一些名称搜索也会返回多个CID。因此，建议永远不要将CID
作为选择列的一部分，因为代码会崩溃。
5. 一些搜索将返回NULL或空列表，原因在结论中有讨论，
但这些问题在准备循环搜索代码时仍需关注。

代码明确处理了上述所有5种情况，如下所述。
"""

new_df = []

for tauro in tauros:
    compounds = pubchem.get_compounds(tauro, "name")
    if compounds:
        # 如果返回多个结果，只保留第一个
        compound = compounds[0]
        compound_dict = compound.to_dict(use_calls)
        
        # 在将其设置为索引之前，检查'cid'是否存在于化合物属性中
        if 'cid' in compound_dict:
            df = pd.DataFrame([compound_dict], index=[compound.cid])
        else:
            df = pd.DataFrame([compound_dict])
        new_df.append(df)

# 合并列表中的所有数据框
tauro_df = pd.concat(new_df, ignore_index=True)
tauro_df['molecular_weight'] = tauro_df['molecular_weight'].astype("float")

# 在数据框中添加mol对象

PandasTools.AddMoleculeColumnToFrame(tauro_df, smilesCol = "canonical_smiles")

#展示信息
tauro_df.info()

#统计信息描述
tauro_df.describe()

# 检查数据分布

nums = ["h_bond_acceptor_count","h_bond_donor_count","molecular_weight","tpsa","xlogp"]
fig, axes = plt.subplots(nrows = 5, ncols = 1, figsize = (9, 15))
plt.tight_layout()
for ax,num in zip(axes,nums):
    ax.hist(tauro_df[num], color = "black")
    ax.set_xlabel(f"{num} in au")
    ax.set_ylabel("frequency")
    
plt.show()

# 现在对完整的列表执行相同的操作

new_df = []

# 获取完整列表中的所有胆汁酸名称
biles = [name for name in acids["RefMet Name"]]

for bile in biles:
    results = pubchem.get_compounds(bile, "name")
    if results:
        # 如果返回多个结果，只保留第一个
        results = results[0]
        result_dict = results.to_dict(use_calls)
        
        # 在将其设置为索引之前，检查'cid'是否存在于化合物属性中
        if 'cid' in result_dict:
            df = pd.DataFrame([result_dict], index=[results.cid])
        else:
            df = pd.DataFrame([result_dict])
        new_df.append(df)

# 合并列表中的所有数据框
final_df = pd.concat(new_df, ignore_index=True)
final_df['molecular_weight'] = final_df['molecular_weight'].astype("float")

# 在数据框中添加mol对象

PandasTools.AddMoleculeColumnToFrame(final_df, smilesCol = "canonical_smiles")

# 输出最终的数据框
final_df

# 检查分布

nums = ["h_bond_acceptor_count", "h_bond_donor_count", "molecular_weight", "tpsa", "xlogp"]
fig, axes = plt.subplots(nrows=5, ncols=1, figsize=(9, 15))
plt.tight_layout()

# 为每个数值属性绘制直方图
for ax, num in zip(axes, nums):
    ax.hist(final_df[num], color="black")
    ax.set_xlabel(f"{num} in au")
    ax.set_ylabel("frequency")
    
plt.show()

代码功能说明：

• 获取完整列表的胆汁酸数据：这段代码首先从acids["RefMet Name"]中获取所有胆汁酸的名称，接着使用pubchem.get_compounds()从PubChem检索相应的化学数据。

• 处理返回的数据：如果返回的结果有多个化合物，只保留第一个，并转换为字典格式。接着，通过判断cid是否存在来设置数据框的索引，避免出现错误。

• 合并数据框：把所有检索到的数据框合并成一个大的数据框final_df，并确保molecular_weight列的数据类型为浮动数值类型。

• 添加分子结构：使用RDKit将每个分子的SMILES字符串转换为分子对象，并将其添加到数据框中。

• 检查属性分布：通过matplotlib绘制多个属性（如氢键受体数、氢键供体数、分子量、拓扑表面积等）的直方图，展示不同分子属性的分布情况。

# 仅在tauro中查找胆甾烷基团；反之，查找区分因素

tauro_mols = [tauro_mol for tauro_mol in tauro_df["ROMol"]]
res_t = rdFMCS.FindMCS(tauro_mols)
res_t.smartsString
# [OUT]'[#6]-[#6](-[#6]-[#6]-[#6](=[#8])-[#7]-[#6]-[#6]-[#16](=[#8])(=[#8])-[#8])-[#6]1-[#6]-[#6]-[#6]2-[#6]-1(-[#6]-[#6]-[#6]1-[#6]-2-[#6]-[#6]-[#6]2-[#6]-1(-[#6]-[#6]-[#6]-[#6]-2)-[#6])-[#6]'

# 在所有胆汁酸中查找胆甾烷基团；反之，查找区分因素

molecules =  [thing for thing in final_df["ROMol"]]
res_m = rdFMCS.FindMCS(molecules)
res_m.smartsString
# [OUT]'[#6]-[#6](-[#6]-[#6])-[#6]1-[#6]-[#6]-[#6]-[#6]-1-[#6]-[#6]-[#6]1-[#6]-[#6]-[#6]-[#6]2-[#6]-1(-[#6]-[#6]-[#6]-[#6]-2)-[#6]'

如果你没有注意到，tauro丰富列表中的识别基团更长。现在，让我们看看这些基团在结构上的样子。

# 显示tauro的基团

res_tmotif = Chem.MolFromSmarts(res_t.smartsString)
res_tmotif

# 显示所有的基团

res_motif = Chem.MolFromSmarts(res_m.smartsString)
res_motif

显著的是，包含牛磺酸和硫酸盐共轭物的尾部主导了tauro基团，而未修改的骨架则主导了综合列表。有趣的是，通常称为C环的第三环在所有胆汁酸的基团中都是开环的，这可能是由于该环中的羟基侧链随机放置所致。尽管如此，这两个代表性基团并不远离预期。你可以在列表中的所有分子中高亮显示该基团（未显示）。

# 在tauro_df中高亮显示基团

highlight_tauro = [tauro_mol.GetSubstructMatch(res_tmotif) for tauro_mol in tauro_mols]
Draw.MolsToGridImage(tauro_mols, highlightAtomLists = highlight_tauro, subImgSize = (250, 250))

你还可以在更大的数组列表中执行：

# 高亮显示所有基团

highlight_mcs = [bile_mol.GetSubstructMatch(res_motif) for bile_mol in molecules]
Draw.MolsToGridImage(molecules, highlightAtomLists = highlight_mcs, subImgSize = (250, 250), useSVG = False, molsPerRow = 4)

输出会有警告，因为显示过多，可自行调整maxMols更改显示数量。

最后，我使用大列表中的13号胆汁酸作为代表，执行了Tanimoto相似性分析。为了节省空间，代码建议的前34个结构的实际输出没有在下方显示。相反，随后的代码片段显示了第13、14和15号结构，它们都是主要的胆汁酸，并通过Tanimoto相似性进行了注释。

from rdkit import DataStructs

# 选择参考分子并执行Tanimoto相似性分析
ref_ECFP4_fps = AllChem.GetMorganFingerprintAsBitVect(final_df['ROMol'][12], 2)

bulk_ECFP4_fps = [AllChem.GetMorganFingerprintAsBitVect(x, 2) for x in final_df['ROMol']]
similarity_efcp4 = [DataStructs.FingerprintSimilarity(ref_ECFP4_fps, x) for x in bulk_ECFP4_fps]

# 我们还可以将similarity_efcp4添加到数据框中，并可视化结构和相似性。
final_df['Tanimoto_Similarity (ECFP4)'] = similarity_efcp4

# 显示前34行
PandasTools.FrameToGridImage(final_df.head(34), legendsCol="Tanimoto_Similarity (ECFP4)", molsPerRow = 4)

# 输出：主要胆汁酸第13号（参考）、第14号和第15号，通过Tanimoto相似性进行了注释。
# 将数据框导出为Excel文件
final_df.to_excel("bile_acids_070224.xlsx")
最后，数据框被导出为Excel电子表格，以便存储和进一步审查（如果需要）。

结果与讨论

根据RDKit中的Mol对象表示，牛磺胆酸共轭胆汁酸在其前后端都经过极性基团（如硫酸盐和氨基酸相关片段）的化学修饰。因此，牛磺胆酸共轭胆汁酸在大小（tpsa）和分子量（MW）上都较大，这在描述性统计和直方图中有所体现。您可以自行验证这一点。
牛磺胆酸共轭胆汁酸还表现出更高的极性，尤其是在氢键供体-受体（HA/HD）度量和亲脂性（XLogP）方面。牛磺胆酸共轭类的亲脂性较低，这在描述性统计和直方图中也有所体现。

根据描述性统计，富集集的总体相似性较高，而整体集的相似性较低。

即使在使用PubChemPy进行精心筛选后，冗余依然可见。需要进一步分析以确定其来源，无论是立体化学问题还是使用了多个冗余标识符。解决此问题超出了本教程的范围，但值得注意。

在使用PubChemPy进行搜索时，列表中的八个标识符返回了空列表。然而，使用Web界面进行手动输入时，返回了以下信息和一些例外（见下）：
手动搜索的例外和错误代码：

• “3,12-二酮胆汁酸” — 错误代码：NULL，

• “3,6-二酮胆汁酸” — 错误代码：“结构不可用”，

• “3,7-二酮胆汁酸” — 该结构是由日本化学物质字典提供的有限审查。尽管如此，暂时可以使用替代CID，直到进一步审查。

• “3-酮-鹅去氧胆酸” — 错误代码：NULL，

• “3-酮-去氧胆酸” — 错误代码：NULL，

• “甘氨牛磺胆酸” — 该结构也来自日本化学物质字典的有限审查。尽管如此，以下标识符可以用来访问该记录，直到进一步审查：物质SID: 274209944，化合物CID: 44263376。

• “牛磺去氧胆酸-3-硫酸盐” — 仅通过Web界面进行手动搜索指向牛磺去氧胆酸-3-硫酸盐-d4（钠盐）。

• “牛磺乌苏去氧胆酸-3-硫酸盐” — 仅通过Web界面进行手动搜索指向牛磺乌苏去氧胆酸-3-硫酸盐钠盐-d4。

本教程展示了使用PubChemPy和RDKit进行基础化学信息学分析的一次练习，分析了胆汁酸等生物重要分子。也许我将本教程定制得如此技术性强，并针对那些在生物化学方面有专门兴趣的人（可以理解为“我自己”）进行，可能使得本练习对于一般读者来说有些不易接近。我表示歉意。不过，我认为这仍然是基础的，因为我想更高级的用户会这样认为，它主要作为一个演示。当我达到那个水平时，或许我会知道如何调整展示方式。现在，对于那些对胆汁酸、RDKit和PubChemPy的技术细节感兴趣的人来说，考虑将其作为一个基础，继续进行扩展。希望有人能像我一样享受这次练习。

原文链接