在化学信息学中，SMILES（Simplified Molecular Input Line Entry System）是一种用于表示分子结构的简洁而人类可读的线性表示法。然而，由于 SMILES 表示法可以有多种不同的形式来表示同一个分子，因此在数据存储和处理时，通常需要一个标准化的表示法。这种标准化的表示法被称为“规范编码”。规范编码确保每个分子都有一个唯一且标准的表示，方便检索和比较。

为了进行标准化分析和数据库存储，通常需要将 SMILES 表示法转换为其对应的规范编码。本篇博客将详细介绍如何使用 Python 库 RDKit 和 NetworkX 来实现这一转换。

一、SMILES 表示法简介

SMILES 是一种用字符串表示化学分子的方式，它将分子的化学结构信息编码为简单的 ASCII 字符串。这种表示法易于阅读和书写，同时也能被计算机程序轻松解析。例如，乙醇的 SMILES 表示法是 CCO，它表示一个包含两个碳（C）和一个氧（O）原子的分子。

二、 规范编码的基本原理

规范编码的基本思想是为分子中的每个原子和键分配一个优先级，然后按这些优先级对分子进行BFS遍历-广度优先遍历算法和 编码。

1 原子和键的优先级定义：

在下图中节点就作为我们原子，而节点与节点之间的符号就是键作为两个节点之间的关系。

原子优先级

precedence_order = {‘P’: 6, ‘B’: 5, ‘I’: 4, ‘S’: 3, ‘N’: 2, ‘O’: 1, ‘C’: 0}

键优先级

bond_symbols = {‘SINGLE’: ‘-’, ‘DOUBLE’: ‘=’, ‘TRIPLE’: ‘≡’, ‘AROMATIC’: ‘:’}
生成分子图：
将分子表示为一个图结构，其中节点代表原子，边代表化学键。这一步通常使用化学信息学库（如 RDKit）来实现。

2. 遍历分子图：

使用广度优先搜索（BFS）遍历算法，首先从优先级最高的原子开始遍历。按照遍历可以得到分支节点，它们之间的先后顺序则通过键优先级排序。

如下图，首先从原中可以发现S的原子优先级级别最高，所以它作为第一个处理的节点。处理的过程我们采用bfs广度优先遍历，遍历的过程便纳入了其他分支节点（分子）这些分支节点排序的方式就按照键的优先级去排序。

最后在遍历过程中记录原子和键的顺序，生成唯一的字符串表示，即规范编码。

详细示例: SMILES 转换规范规则

转换步骤和规则

1. 初始节点选择：从所有节点中选择优先级最高的节点作为初始节点。
2. 子节点选择：计算每个相邻节点的优先级和边优先级，选择优先级最高的作为下一个节点。
3. 处理相同优先级的节点和边：如果相邻节点和边的优先级相同，保留所有可能的分支树，直到找到不同的节点。

具体示例

假设我们有以下 SMILES 表示：

O=[N+]([O-])[O-]

初始节点选择

N 是优先级最高的节点，因此选择 N 作为初始节点。

Canonical: N

子节点选择

N 相邻的节点有 O 和 O-，边类型都是单键。

Canonical: N 0-O1 0-O2 0=O3

相同优先级处理

如果相邻节点和边的优先级相同，保留所有可能的分支树，直到找到不同的节点。例如，在选择 O 的子节点时，如果有多个 O，继续搜索直到找到不同的优先级节点。

生成规范表示

最终生成的规范表示如下：

Canonical: N 0-O1 0-O2 0=O3

安装依赖

Conda 环境中，安装以下主要依赖包：

conda create -n GBML pyg=2.4.0 python=3.8 matplotlib python rdkit networkx jupyter pytorch torchvision torchaudio cpuonly pytorch-scatter pytorch-sparse -c pytorch -c pyg -c conda-forge -y

1. SMILES 到 NetworkX 图的转换

我们将使用 RDKit 将 SMILES 表示法转换为分子对象，并使用 NetworkX 将其表示为图结构。在这个图中，节点表示原子，边表示化学键。

import networkx as nx
from rdkit import Chem

def smiles_to_nx_mol(smiles):
    mol = Chem.MolFromSmiles(smiles)
    Chem.Kekulize(mol)
    bond_symbols = {'SINGLE': '-', 'DOUBLE': '=', 'TRIPLE': '≡', 'AROMATIC': ':'}

    G = nx.Graph()

    for atom in mol.GetAtoms():
        G.add_node(atom.GetIdx(), atom_symbol=atom.GetSymbol())

    for bond in mol.GetBonds():
        G.add_edge(bond.GetBeginAtomIdx(),
                   bond.GetEndAtomIdx(),
                   bond_type=bond_symbols[str(bond.GetBondType()).split('.')[-1]])

    return G, mol

2. 定义原子和键的优先级

为了标准化表示，我们需要定义原子和键的优先级。

def get_atom_precedence(G):
    precedence_order = {'P': 6, 'B': 5, 'I': 4, 'S': 3, 'N': 2, 'O': 1, 'C': 0}
    precedence = {atom: precedence_order.get(atom, len(precedence_order) + 1) for atom in
                  nx.get_node_attributes(G, 'atom_symbol').values()}
    return precedence

def get_bond_precedence(bond_type):
    bond_order = {'-': 1, '=': 2, '≡': 3, ':': 4}
    return bond_order.get(bond_type, 0)

3.广度优先搜索（BFS）遍历

使用广度优先搜索遍历图结构，生成规范编码。我们使用 BFS 算法遍历图。在每一步中，按照节点和键的优先级选择下一个节点。如果多个节点和键的优先级相同，则保留所有可能的分支，直到找到唯一的路径。

from collections import deque

def bfs_traversal(start_node, sorted_nodes, G, atom_precedence, bond_priority):
    visited = set()
    node_order = {}
    edge_order = []
    queue = deque([(start_node, None)])
    atom_index = 0
    previous_node = None

    while queue or sorted_nodes:
        if not queue and sorted_nodes:
            while sorted_nodes:
                next_node = sorted_nodes.pop(0)
                if next_node not in visited:
                    queue.append((next_node, previous_node))  # 回到上一个节点的上下文
                    break
            if not queue:
                break

        current_node, parent = queue.popleft()
        if current_node not in visited:
            visited.add(current_node)
            node_order[current_node] = atom_index
            atom_index += 1

            neighbors = list(G.neighbors(current_node))
            neighbors = [n for n in neighbors if n not in visited]

            # 先按原子优先级排序
            neighbors.sort(key=lambda x: -atom_precedence[G.nodes[x]['atom_symbol']])
            # 再按键的优先级和原子优先级综合排序
            neighbors.sort(key=lambda x: (bond_priority[G.edges[current_node, x]['bond_type']],
                                          -atom_precedence[G.nodes[x]['atom_symbol']], x))

            # 检查是否所有邻居的原子和键优先级都相同
            if len(neighbors) > 1 and all(
                    atom_precedence[G.nodes[neighbors[0]]['atom_symbol']] == atom_precedence[G.nodes[n]['atom_symbol']] and
                    bond_priority[G.edges[current_node, neighbors[0]]['bond_type']] == bond_priority[G.edges[current_node, n]['bond_type']]
                    for n in neighbors[1:]):
                # 更深层次的比较，根据邻居的邻居的原子和键优先级进行进一步排序
                neighbors.sort(key=lambda x: (
                    sorted((atom_precedence[G.nodes[n]['atom_symbol']], bond_priority[G.edges[x, n]['bond_type']])
                           for n in G.neighbors(x) if n not in visited),
                    x
                ))

            for neighbor in neighbors:
                if neighbor not in visited:
                    queue.append((neighbor, current_node))
                    bond = G.edges[current_node, neighbor]['bond_type']
                    edge_order.append((node_order[current_node], bond, neighbor))  # 保持邻居索引不变

            previous_node = current_node  # 更新上一个节点

    return node_order, edge_order

4.生成规范编码

根据 BFS 遍历的结果，生成最终的规范表示。这包括原子的顺序和键的类型，确保唯一性和标准化。定义一个函数来整合以上步骤，生成最终的规范编码。

def find_canonical(smiles):
    G, mol = smiles_to_nx_mol(smiles)

    atom_precedence = get_atom_precedence(G)
    bond_priority = {'-': 1, '=': 2, '≡': 3, ':': 4}

    nodes = list(G.nodes(data=True))
    nodes.sort(key=lambda x: (-atom_precedence[x[1]['atom_symbol']], x[0]))

    highest_priority = -atom_precedence[nodes[0][1]['atom_symbol']]
    highest_priority_nodes = [node for node in nodes if -atom_precedence[node[1]['atom_symbol']] == highest_priority]

    if len(highest_priority_nodes) > 1:
        highest_priority_nodes.sort(key=lambda x: -len(list(G.neighbors(x[0]))))
        max_neighbors = len(list(G.neighbors(highest_priority_nodes[0][0])))
        same_neighbors_nodes = [node for node in highest_priority_nodes if len(list(G.neighbors(node[0]))) == max_neighbors]

        if len(same_neighbors_nodes) > 1:
            same_neighbors_nodes.sort(key=lambda x: (
                max(bond_priority[G.edges[x[0], neighbor]['bond_type']] for neighbor in G.neighbors(x[0])),
                x[0]
            ))
            highest_priority_nodes = same_neighbors_nodes

    remaining_nodes = [node for node in nodes if node not in highest_priority_nodes]
    nodes = highest_priority_nodes + remaining_nodes

    best_start_node = nodes[0][0]
    node_order, edge_order = bfs_traversal(best_start_node, [n[0] for n in nodes], G, atom_precedence, bond_priority)

    final_code_parts = [f"{G.nodes[best_start_node]['atom_symbol']}"]
    for parent, bond, child in edge_order:
        final_code_parts.append(f"{parent}{bond}{G.nodes[child]['atom_symbol']}{node_order[child] }")

    final_code = " ".join(final_code_parts)

    return final_code

5. 测试函数

使用一些测试用例来验证我们的函数是否工作正常.

if __name__ == "__main__":
    import pandas as pd

    test_cases = pd.read_csv("test_data_task1.csv")

    for idx, row in test_cases.iterrows():
        smiles = row['smiles']
        expected_output = row['canonical']
        actual_output = find_canonical(smiles)

        print(f"\n{idx + 1} 测试 SMILES: {smiles}")
        print(f"期望输出: {expected_output}")
        print(f"实际输出: {actual_output}")

        if actual_output == expected_output:
            print("✅")
        else:
            print("❌")

3.整体代码

from collections import deque
import networkx as nx
from rdkit import Chem

def smiles_to_nx_mol(smiles):
    mol = Chem.MolFromSmiles(smiles)
    Chem.Kekulize(mol)
    bond_symbols = {'SINGLE': '-', 'DOUBLE': '=', 'TRIPLE': '≡', 'AROMATIC': ':'}

    G = nx.Graph()

    for atom in mol.GetAtoms():
        G.add_node(atom.GetIdx(), atom_symbol=atom.GetSymbol())

    for bond in mol.GetBonds():
        G.add_edge(bond.GetBeginAtomIdx(),
                   bond.GetEndAtomIdx(),
                   bond_type=bond_symbols[str(bond.GetBondType()).split('.')[-1]])

    return G, mol

def get_atom_precedence(G):
    precedence_order = {'P': 6, 'B': 5, 'I': 4, 'S': 3, 'N': 2, 'O': 1, 'C': 0}
    precedence = {atom: precedence_order.get(atom, len(precedence_order) + 1) for atom in
                  nx.get_node_attributes(G, 'atom_symbol').values()}
    return precedence

def get_bond_precedence(bond_type):
    bond_order = {'-': 1, '=': 2, '≡': 3, ':': 4}
    return bond_order.get(bond_type, 0)

def bfs_traversal(start_node, sorted_nodes, G, atom_precedence, bond_priority):
    visited = set()
    node_order = {}
    edge_order = []
    queue = deque([(start_node, None)])
    atom_index = 0
    previous_node = None

    while queue or sorted_nodes:
        if not queue and sorted_nodes:
            while sorted_nodes:
                next_node = sorted_nodes.pop(0)
                if next_node not in visited:
                    queue.append((next_node, previous_node))  # 回到上一个节点的上下文
                    break
            if not queue:
                break

        current_node, parent = queue.popleft()
        if current_node not in visited:
            visited.add(current_node)
            node_order[current_node] = atom_index
            atom_index += 1

            neighbors = list(G.neighbors(current_node))
            neighbors = [n for n in neighbors if n not in visited]

            # 先按原子优先级排序
            neighbors.sort(key=lambda x: -atom_precedence[G.nodes[x]['atom_symbol']])
            # 再按键的优先级和原子优先级综合排序
            neighbors.sort(key=lambda x: (bond_priority[G.edges[current_node, x]['bond_type']],
                                          -atom_precedence[G.nodes[x]['atom_symbol']], x))

            # 检查是否所有邻居的原子和键优先级都相同
            if len(neighbors) > 1 and all(
                    atom_precedence[G.nodes[neighbors[0]]['atom_symbol']] == atom_precedence[G.nodes[n]['atom_symbol']] and
                    bond_priority[G.edges[current_node, neighbors[0]]['bond_type']] == bond_priority[G.edges[current_node, n]['bond_type']]
                    for n in neighbors[1:]):
                # 更深层次的比较，根据邻居的邻居的原子和键优先级进行进一步排序
                neighbors.sort(key=lambda x: (
                    sorted((atom_precedence[G.nodes[n]['atom_symbol']], bond_priority[G.edges[x, n]['bond_type']])
                           for n in G.neighbors(x) if n not in visited),
                    x
                ))

            for neighbor in neighbors:
                if neighbor not in visited:
                    queue.append((neighbor, current_node))
                    bond = G.edges[current_node, neighbor]['bond_type']
                    edge_order.append((node_order[current_node], bond, neighbor))  # 保持邻居索引不变


            previous_node = current_node  # 更新上一个节点


    return node_order, edge_order

def find_canonical(smiles):
    G, mol = smiles_to_nx_mol(smiles)

    atom_precedence = get_atom_precedence(G)
    bond_priority = {'-': 1, '=': 2, '≡': 3, ':': 4}

    nodes = list(G.nodes(data=True))


    edges = list(G.edges(data=True))


    # 首次排序，根据原子优先级
    nodes.sort(key=lambda x: (-atom_precedence[x[1]['atom_symbol']], x[0]))


    # 检查是否有多个原子具有相同的最高优先级
    highest_priority = -atom_precedence[nodes[0][1]['atom_symbol']]
    highest_priority_nodes = [node for node in nodes if -atom_precedence[node[1]['atom_symbol']] == highest_priority]

    if len(highest_priority_nodes) > 1:
        # 如果有多个最高优先级的原子，则根据邻居数量进行排序（降序）
        highest_priority_nodes.sort(key=lambda x: -len(list(G.neighbors(x[0]))))

        # 检查是否有多个节点具有相同的邻居数量
        max_neighbors = len(list(G.neighbors(highest_priority_nodes[0][0])))
        same_neighbors_nodes = [node for node in highest_priority_nodes if len(list(G.neighbors(node[0]))) == max_neighbors]

        if len(same_neighbors_nodes) > 1:
            # 如果有多个节点具有相同的邻居数量，则根据键的优先级排序
            same_neighbors_nodes.sort(key=lambda x: (
                max(bond_priority[G.edges[x[0], neighbor]['bond_type']] for neighbor in G.neighbors(x[0])),
                x[0]
            ))
            highest_priority_nodes = same_neighbors_nodes



    # 将最高优先级节点列表与剩余节点列表合并
    remaining_nodes = [node for node in nodes if node not in highest_priority_nodes]
    nodes = highest_priority_nodes + remaining_nodes

    # 打印重新排序后的完整节点列表

    # 选择优先级最高的节点作为起始节点
    best_start_node = nodes[0][0]


    # 开始 BFS 遍历
    node_order, edge_order = bfs_traversal(best_start_node, [n[0] for n in nodes], G, atom_precedence, bond_priority)

    # 生成最终的编码
    final_code_parts = [f"{G.nodes[best_start_node]['atom_symbol']}"]
    for parent, bond, child in edge_order:
        final_code_parts.append(f"{parent}{bond}{G.nodes[child]['atom_symbol']}{node_order[child] }")

    final_code = " ".join(final_code_parts)


    return final_code

if __name__ == "__main__":
    import pandas as pd

    test_cases = pd.read_csv("test_data_task1.csv")

    for idx, row in test_cases.iterrows():
        smiles = row['smiles']
        expected_output = row['canonical']
        actual_output = find_canonical(smiles)

        print(f"\n{idx + 1} 测试 SMILES: {smiles}")
        print(f"期望输出: {expected_output}")
        print(f"实际输出: {actual_output}")

        if actual_output == expected_output:
            print("✅")
        else:
            print("❌")

导入测试文件.csv

smiles,canonical
O=[N+]([O-])[O-],N 0-O1 0-O2 0=O3
O=S(=O)([O-])[O-],S 0-O1 0-O2 0=O3 0=O4
O=S([O-])([O-])=S,S 0-O1 0-O2 0=S3 0=O4
C#CC1(OC(N)=O)CCCCC1,N 0-C1 1-O2 1=O3 2-C4 4-C5 4-C6 4-C7 5-C8 6-C9 7≡C10 8-C11 9-C11
O=C([O-])C(S)C(S)C(=O)[O-],S 0-C1 1-C2 1-C3 2-S4 2-C5 3-O6 3=O7 5-O8 5=O9
O=C([O-])/C=C/C(=O)[O-],O 0-C1 1-C2 1=O3 2=C4 4-C5 5-O6 5=O7
O=C([O-])CCC(=O)[O-],O 0-C1 1-C2 1=O3 2-C4 4-C5 5-O6 5=O7
O=C([O-])CCCCCCCC(=O)[O-],O 0-C1 1-C2 1=O3 2-C4 4-C5 5-C6 6-C7 7-C8 8-C9 9-C10 10-O11 10=O12
O=S(=O)([O-])CCS,S 0-O1 0-C2 0=O3 0=O4 2-C5 5-S6
O=C([O-])P(=O)([O-])[O-],P 0-O1 0-O2 0-C3 0=O4 3-O5 3=O6

五、总结
通过上述步骤，我们可以将 SMILES 表示法转换为规范编码。这种标准化的表示法不仅有助于数据存储和检索，还可以用于分子相似性比较和其他化学信息学分析。

4. 测试结果

1 测试 SMILES: O=N+[O-]
期望输出: N 0-O1 0-O2 0=O3
实际输出: N 0-O1 0-O2 0=O3
✅

2 测试 SMILES: O=S(=O)([O-])[O-]
期望输出: S 0-O1 0-O2 0=O3 0=O4
实际输出: S 0-O1 0-O2 0=O3 0=O4
✅

3 测试 SMILES: O=S([O-])([O-])=S
期望输出: S 0-O1 0-O2 0=S3 0=O4
实际输出: S 0-O1 0-O2 0=S3 0=O4
✅

4 测试 SMILES: C#CC1(OC(N)=O)CCCCC1
期望输出: N 0-C1 1-O2 1=O3 2-C4 4-C5 4-C6 4-C7 5-C8 6-C9 7≡C10 8-C11 9-C11
实际输出: N 0-C1 1-O2 1=O3 2-C4 4-C5 4-C6 4-C7 5-C8 6-C9 7≡C10 8-C11 9-C11
✅

5 测试 SMILES: O=C([O-])C(S)C(S)C(=O)[O-]
期望输出: S 0-C1 1-C2 1-C3 2-S4 2-C5 3-O6 3=O7 5-O8 5=O9
实际输出: S 0-C1 1-C2 1-C3 2-S4 2-C5 3-O6 3=O7 5-O8 5=O9
✅

6 测试 SMILES: O=C([O-])/C=C/C(=O)[O-]
期望输出: O 0-C1 1-C2 1=O3 2=C4 4-C5 5-O6 5=O7
实际输出: O 0-C1 1-C2 1=O3 2=C4 4-C5 5-O6 5=O7
✅

7 测试 SMILES: O=C([O-])CCC(=O)[O-]
期望输出: O 0-C1 1-C2 1=O3 2-C4 4-C5 5-O6 5=O7
实际输出: O 0-C1 1-C2 1=O3 2-C4 4-C5 5-O6 5=O7
✅

8 测试 SMILES: O=C([O-])CCCCCCCC(=O)[O-]
期望输出: O 0-C1 1-C2 1=O3 2-C4 4-C5 5-C6 6-C7 7-C8 8-C9 9-C10 10-O11 10=O12
实际输出: O 0-C1 1-C2 1=O3 2-C4 4-C5 5-C6 6-C7 7-C8 8-C9 9-C10 10-O11 10=O12
✅

9 测试 SMILES: O=S(=O)([O-])CCS
期望输出: S 0-O1 0-C2 0=O3 0=O4 2-C5 5-S6
实际输出: S 0-O1 0-C2 0=O3 0=O4 2-C5 5-S6
✅

10 测试 SMILES: O=C([O-])P(=O)([O-])[O-]
期望输出: P 0-O1 0-O2 0-C3 0=O4 3-O5 3=O6
实际输出: P 0-O1 0-O2 0-C3 0=O4 3-O5 3=O6
✅