如何使用RDkit计算Morgan指纹并对其进行可视化

Morgan指纹中更广为人知的是圆形指纹，该指纹是通过将 Morgan 算法应用于一组用户提供的原子不变量而构建的（具体原理参考上一条帖子）。 生成 Morgan 指纹时，还必须提供指纹的半径。Morgan指纹，例如原子对（atom pair）和拓扑扭转（topological torsions），默认使用计数，但也可以将它们计算为位向量。下面我们以丁酰胺为例向大家展示如何计算其指纹（计数和位相连）：

图1 丁酰胺的结构式

from rdkit.Chem import AllChem
from rdkit import Chem
m1 = Chem.MolFromSmiles('CCCC(N)=O')
fp1 = AllChem.GetMorganFingerprint(m1,2,useCounts=True)
# <rdkit.DataStructs.cDataStructs.UIntSparseIntVect object at 0x...>
# 得到一个需要转换的 UIntSparseIntVect（整数向量）
fp2= AllChem.GetMorganFingerprintAsBitVect(m1,2,nBits=1024)
#<rdkit.DataStructs.cDataStructs.ExplicitBitVect object at 0x...>
# 得到一个需要转换的ExplicitBitVect（1024位向量）
fp3= AllChem.GetMorganFingerprintAsBitVect(m1,2)
#<rdkit.DataStructs.cDataStructs.ExplicitBitVect object at 0x...>
# 得到一个需要转换的ExplicitBitVect（2048位向量）
fp2.ToBitString() #00000000，将fp2转换为位串
fp3.ToBitString()  #00000000，将fp3转换为位串

图2 丁酰胺的特征在去除重复后所代表的结构。该分子在初始赋值阶段产生五个特征（即唯一的原子标识符），在一次迭代后产生六个额外的特征，在两次迭代后产生三个特征。进一步的迭代没有产生新的特征。(图中的中心原子通过显示其原子编号来表示）

关于ECFP4的计数表示方法，我们参考文献1中的描述，当丁酰胺经过两次迭代（n=2），即可以得到ECFP4的3个特征，

图3 不同直径的丁酰胺指纹。 请注意，较大直径包含所有较小直径的指纹位，可能在末尾附加了新标识符。 另外，请注意 ECFP_4 和 ECFP_6 包含相同的列表。 这是因为最后的迭代没有发现任何新的标识符，这里的 “新”是由构成特定特征的键的集合决定的。 当我们达到最大直径为四个键的时候，整个分子已被覆盖，没有什么新的东西可以发现了。

然后我们尝试用代码复现上述ECFP4的计数表示方法，得出丁酰胺的SparseIntVects

print(fp1.GetLength())  # 查看长度，返回指纹的长度
# 4294967295
len(fp1.GetNonzeroElements()) # 查看非空元素的长度
# 14
print(fp1.GetNonzeroElements()) 
# 查看非空元素，返回值是一个字典，键对应非空指纹的位数，值是指纹出现的频数
# {474023245:1, 847957139:1, 864942730:1, 1167322652:1, 1173125914:1, 1510328189:1, 2231929377:1, 2245384272:2, 2246699815:1, 2246728737:1, 2591432844:1, 3328145258:1, 3542456614:1, 4112692162:1}
# 我们可以看到共有14个非空指纹位数，数量与文献中（图3，ECFP4）报道一致，但是具体数值不一样，我猜测可能是编码的方式不一样，但是其代表的结构应该是一样的。
# 接下来我们看看每个非空指纹位数代表的什么。
info={}
fp1 = AllChem.GetMorganFingerprint(m1,2,bitInfo=info)
len(info)
# 14
info[474023245]
# ((1,2), )
#位 474023245由原子1设置了一次，半径为2。
info[2245384272]
# ((1,0), (2,0))
位 2245384272 设置了两次，一次由原子1，一次由原子2，半径为0。
info[4112692162]
# ((2,2), )
# 位4112692162由原子2设置了一次，半径为2。

接下来我们关注474023245位，我们可以提取出原子1的半径为2内的所有原子组成的亚分子

env = Chem.FindAtomEnvironmentOfRadiusN(m1,2,1)
amap={}
submol=Chem.PathToSubmol(m1,env,atomMap=amap)
submol.GetNumAtoms()
#得到原子数量，结果为 4
amap
# {0: 0, 1: 1, 2: 2, 3: 3}，我查了很多资料，还是不明白这个字典代表的意思，只知道键代表的是原子编号
Chem.MolToSmiles(submol)
# 我们通过为该“子分子”生成 SMILES 来“解释”该位
# CCCC

图4 我们将这14个子结构与图2中的理论结构进行一一对应，红色等号的是使用以上代码生成的结构，红色等号右边的是图2种的理论结构。

以上展现了如何根据count指纹生成图片，那么如何根据位指纹生成图片呢？

from rdkit.Chem import Draw
mol = Chem.MolFromSmiles('CCCC(N)=O')  #还是以丁酰胺为例子
bi = {}
fp = AllChem.GetMorganFingerprintAsBitVect(mol, radius=2, bitInfo=bi)
print(bi)  
# 生成一个包含14个键的字典，这14个键分别代表图2和图4中的14个片段
# {80: ((1, 0), (2, 0)), 140: ((4, 1),), 294: ((0, 1),), 545: ((3, 2),), 650: ((5, 0),), 794: ((1, 1),), 807: ((3, 0),), 962: ((2, 2),), 1057: ((0, 0),), 1171: ((4, 0),), 1357: ((1, 2),), 1564: ((2, 1),), 1898: ((3, 1),), 1917: ((5, 1),)}
# 以上键代表2048位中的位名
bi[140] 
#假设我们想显示键140对应的结构
# ((4, 1),)，结果显示为以原子4为中心，以1为半径
Draw.DrawMorganBit(mol, 140, bi, useSVG=True)
# 如下图所示，与前面基于count的指纹能对的上

参考：

1. Rogers, David Hahn, Mathew. Extended-connectivity fingerprints. Journal of Chemical Information and Modeling. 2010.50(5).742-754.
2. https://www.rdkit.org/docs/GettingStartedInPython.html