RDkit一：利用RDkit筛选化学小分子结构特征

        最近在做关于小分子筛选的项目，涉及到一些处理方式，后面会逐篇总结遇到的一些问题及解决方法。

        首先简单介绍一下RDkit，RDkit是一款用来处理化学小分子的python开源包，最早是由诺华公司开发的，底层是由C++编写的。在github上有它的源码和说明书，地址如下：

https://github.com/rdkit/rdkithttps://github.com/rdkit/rdkit        RDkit在Anaconda或miniconda环境中的安装：

conda install -yq -c rdkit rdkit

        再来简单介绍一下smiles编码，smiles编码其实本质上就是一种用字符串表达小分子二维结构的编码，通过rdkit包，可以将chemdraw中.mol格式的小分子结构式，或sdf格式的空间坐标转化成smiles编码，从而带入到机器学习和深度学习模型中进行学习。

下面要进入正题了：

        在一百万个化学小分子中，进行筛选，仅留下NO2<2个，Cl < 3个，Br < 2个，F < 6个，芳环数量 < 5个的分子。

为了解决这个问题，需要用到两个RDkit包，这里只介绍最简单的用法，感兴趣的可以去看包简介或源码：

1. rdkit.Chem.Lipinski

Lipinski(利平普斯规则)是一种小分子药物的常见约束规则，在rdkit的Lipinski包中包含多种参数的计算，只要知道小分子的smiles编码，就可以解析出其，HeavyAtomCount，NumAromaticRings……等多种信息，其中比较常用的包括:NumAromaticRings,NumHAcceptors,NumHDonors,NumRotatableBonds这几种。

2.rdkit.Chem

rdkit.Chem包中含有操作小分子对象的功能，包含原子操作，键操作，环操作，药效团查找等这些功能。

我们这里主要进行的原子操作，原子操作涉及到函数包括：

• 对原子进行遍历：m.GetAtoms()
• 获取原子索引：GetIdx()
• 获取原子序号：GetAtomicNum()
• 获取原子符号：GetSymbol()
• 获取原子连接数（受H是否隐藏影响）：GetDegree()
• 获取原子总连接数（与H是否隐藏无关）：GetTotalDegree()
• 获取原子形式电荷：GetFormalCharge()
• 获取原子杂化方式：GetHybridization()
• 获取原子显式化合价：GetExplicitValence()
• 获取原子隐式化合价：GetImplicitValence()
• 获取原子总的化合价：GetTotalValence()

有一篇知乎文章可以看，写的很好：

RDKit|分子基础操作与药效团查找 - 知乎文章目录 原子操作键操作环操作手动实现氧族药效团查找1.原子操作在rdkit中，分子中的每一个原子都是对象，可以通过原子对象的属性和函数来获取各种信息。 对原子进行遍历：m.GetAtoms()获取原子索引：GetIdx()获…https://zhuanlan.zhihu.com/p/143111689在开始上代码之前，再介绍一下输入数据的格式，因为我发现对于一些接触RDkit比较晚的老生化人来说，还是要写的细致一些才好。输入数据是分子的smiles编码，如下，将所有需要筛选分子的smiles编码放在一个csv或txt文件中即可，列名为smiles：

 不啰嗦，上代码：

import pandas as pd
from rdkit import Chem
from rdkit.Chem import Lipinski

# 读入数据
df = pd.read_csv('smiles.csv')

# 筛选芳环数目：aromatic rings < 5
df['NumAromaticRings'] = df['smiles'].apply(lambda x:Lipinski.NumAromaticRings(Chem.MolFromSmiles(x)))
# 也可并行，用parallel_apply，根据自己需要
df = df.drop(df[df.NumAromaticRings >= 5].index)

# 筛选F,Cl,Br等元素的信息
m = [Chem.MolFromSmiles(x) for x in df.smiles.tolist()]
# 得到氟数目列表
num_F = []
for i in range(len(m)):
    F = [atom.GetSymbol() for atom in m[i].GetAtoms()].count('F')
    num_F.append(F)

# 得到含氯数目列表
num_Cl = []
for i in range(len(m)):
    Cl = [atom.GetSymbol() for atom in m[i].GetAtoms()].count('Cl')
    num_Cl.append(Cl)

# 得到含溴数目列表
num_Br = []
for i in range(len(m)):
    Br = [atom.GetSymbol() for atom in m[i].GetAtoms()].count('Br')
    num_Br.append(Br)

# F,Br,CL数量均可筛选
print(f'max_F:{max(num_F)};max_Cl:{max(num_Cl)};max_Br:{max(num_Br)}')

# 条件筛选
transform = {'num_F':num_F,'num_Cl':num_Cl,'num_Br':num_Br}
atom_nums = pd.DataFrame(transform)
df = pd.concat([df,atom_nums],axis=1)
print(df.info())
df = df.drop(df[df.num_F >= 6].index)
df = df.drop(df[df.num_Cl >= 3].index)
df = df.drop(df[df.num_Br >= 2].index)

# 筛选NO2基团，此处是根据字符串串查找进行的
dff = df[df['smiles'].str.contains(pat='[N+](=O)[O-]',regex=False)]
# 得到含NO2基团的smiles编码dataframe表格
smiles = dff.smiles.tolist()
NO2 = [xcount('[N+](=O)[O-]') for x in n]    
transform2 = {'smiles':smiles,'num_NO2':num_NO2}
FG_nums = pd.DataFrame(transform2)
# 以smiles为inner与df求差集
df2 = df.append(FG_nums)
df2 = df2.drop_duplicates(subset=['smiles'],keep=False)

# 保存最终结果
df2.to_csv('result.csv',header = True, index = False)

每周一篇，下周预告RDkit处理聚类问题