RDKit|化学特征、药效团提取与2D药效团指纹计算

  • 一、化学特征和药效团提取
  • 二、化学特征文件介绍
    • 1.化学特征(chemical features)
    • 2.FDef文件语法
  • 三、2D药效团指纹
    • 1.编码原理
    • 2.参数设置
    • 3.生成2D药效团指纹
    • 4.修改FDef设置
    • 5.Gobbi 2D药效团指纹

一、化学特征和药效团提取

Rdkit中有一系列基于SMARTS定义的化学特征。要对分子进行化学特征分析和计算时,需要先导入一个特征库,创建一个特征工厂,并通过特征工厂计算化学特征。先来操作一下,原理及介绍在第二部分。

  • 获取特征库:RDConfig.RDDataDir目录下的'BaseFeatures.fdef'
  • 构建特征工厂:ChemicalFeatures.BuildFeatureFactory(fdefName)
    fdefName:特征库文件
>>> from rdkit import Chem
>>> from rdkit.Chem import ChemicalFeatures
>>> from rdkit import RDConfig
>>> import os
>>> fdefName = os.path.join(RDConfig.RDDataDir,'BaseFeatures.fdef')
>>> factory = ChemicalFeatures.BuildFeatureFactory(fdefName)
  • 使用特征工厂搜索特征:GetFeaturesForMol(m)
>>> m = Chem.MolFromSmiles('OCc1ccccc1CN')
>>> feats = factory.GetFeaturesForMol(m)
>>> print(len(feats))
8

搜索到的每个特征都包含了该特征家族(例如供体、受体等)、特征类别、该特征对应的原子、特征对应序号等信息。

  • 特征家族信息:GetFamily()
  • 特征类型信息:GetType()
  • 特征对应原子:GetAtomIds()
  • 特征对应序号:GetId()
>>> print(feats[0].GetFamily())
>>> print(feats[0].GetType())
>>> print(feats[0].GetAtomIds())
Donor
SingleAtomDonor
(0,)
  • 如果分子包含坐标信息,化学特征也会包括原子坐标:GetPos()
>>> from rdkit.Chem import AllChem
>>> AllChem.Compute2DCoords(m)
>>> list(feats[0].GetPos())
[2.0705367611607857, -2.335674960409047, 0.0]

二、化学特征文件介绍

以上都是基于现有的特征库进行分析和提取,而特征库就是一个特征定义文件(Feature Definition File, FDef)。该文件包含了一系列的化学特征及它们的所有信息,并通过SMARTS来表示。除了化学特征,FDef文件也有对原子类型的定义及解释,让特征更容易理解。

  • 每个化学特征对应的SMARTS,对SMARTS不了解的可以参考这份速查表
    key为“特征家族.特征类型”,值为SMARTS
>>> list(factory.GetFeatureDefs().items())[:2]
[('Donor.SingleAtomDonor',
  '[$([N&!H0&v3,N&!H0&+1&v4,n&H1&+0,$([$([Nv3](-C)(-C)-C)]),$([$(n[n;H1]),$(nc[n;H1])])]),$([O,S;H1;+0])]'),
 ('Acceptor.SingleAtomAcceptor',
  '[$([O;H1;v2]),$([O;H0;v2;!$(O=N-*),$([O;-;!$(*-N=O)]),$([o;+0])]),$([n;+0;!X3;!$([n;H1](cc)cc),$([$([N;H0]#[C&v4])]),$([N&v3;H0;$(Nc)])]),$([F;$(F-[#6]);!$(FC[F,Cl,Br,I])])]')]

1.化学特征(chemical features)

化学特征由一个特征类型(Feature Type)和特征家族(Feature Family)共同定义。特征家族是对特征总体上的分类,例如氢键供体、芳香性等。药效团匹配就是根据特征家族来实现的。而特征类型提供了关于特征的一些更详细的信息。每个特征类型包含了以下信息:

  • 一个SMARTS表达式,定义了该特征类型所匹配的原子。
  • 特征位置的权重,特征位置由原子位置决定。

2.FDef文件语法

FDef中包含两个概念:原子类型和特征。可以大致这么认为:原子类型是最底层的抽象,对相似的原子做了归类。特征是更高级的抽象,对相似的原子和原子类型再进行归类。以下都围绕原子类型定义和特征定义的语法展开。原文在这里

(1)原子类型(AtomType)定义

  • 原子类型定义相当于给用于查询原子的SMARTS起了一个别名,可以让FDef更具可读性。例如,可以用如下方式定义一个非极性碳原子,Carbon_NonPolar就是一个原子类型定义名称:
    AtomType Carbon_NonPolar [C&!$(C=[O,N,P,S])]

  • 因此,可以创建一个有利于理解和使用SMARTS的名称。要引用某个原子类型,可以写下该名称并用大括号括起来。例如,定义另一个原子类型,取名为Hphobe,让Hphobe包含Carbon_NonPolar:
    AtomType Carbon_NonPolar [C&!$(C=[O,N,P,S])]
    AtomType Hphobe [{Carbon_NonPolar},c,s,S&H0&v2,F,Cl,Br,I]

  • 重复写下一个原子类型时,意味着把他们合在一起,相当于用“,”连接(在SMARTS中表示或),例如下面两个语句:
    AtomType d1 [N&!H0]
    AtomType d1 [O&!H0]

  • 等价于:
    AtomType d1 [N&!H0,O&!H0]

  • 更简洁的写法:
    AtomType d1 [N,O;!H0]

  • 要注意“&”和“;”都表示“与”,“,”表示“或”,但“&”的优先级高于“,”,而“;”的优先级低于“,”。

  • 类似与SMARTS,原子特征定义也可以用“!”来表示非,而“!”会与自身的表达式结合,例如:
    AtomType d1 [N,O,S]
    AtomType !d1 [H0]

  • “!d1”等价于:
    AtomType d1 [!H0;N,O,S]

(2)特征(Feature)定义

  • 特征定义比原子类型定义更复杂,并且需要多行实现,例如:
    DefineFeature HDonor1 [N,O;!H0]
    Family HBondDonor
    Weights 1.0
    EndFeature

  • 特征定义的第一行包含了特征类型和所规定的SMARTS,第二行和第三行(没有先后顺序)定义了特征家族和原子权重(权重值的数量和SMARTS中包含的原子数量相同,有多个权重时,用逗号分隔)。原子权重用来计算特征的位置。最后一行表示结束定义,必须是“EndFeature”。原子类型定义和特征定义可以混在一起使用,只要在引用前定义好就行了。

(3)其它语法

  • 井号“#”开头表示注释,程序运行时会忽略该行。
  • 反斜杠“\”结尾表示续行,也就是该行没有结束。
  • 开头的空格将被忽略。

三、2D药效团指纹

1.编码原理

对2D药效团的解释如下图。这里做一个简单的解释。2D药效团指纹,是在药效团的基础上,加上2D(拓扑)距离实现。当在药效团间加入了距离信息后,每个特征(药效团+距离)会被分配一个唯一的id,所有的id共同组成了分子的2D药效团指纹。指纹生成过程中,有两个关键参数,第一个是组成指纹的药效团数量,需要提供最少和最多药效团数。下图中的“2-3 point pharmacophores”表示一个指纹特征中,最少包括2个,最多包括3个药效团。第二个关键参数是拓扑距离区间,用于对两个药效团间的距离进行归类,对应下图的“distance bins”。根据药效团组成和距离的不同,形成不同的bit id。在Rdkit的体系中,药效团指纹的bit ids是没有经过哈希映射的,也就是说,指纹的每一位都代表了一组特征和它们的距离。所以,给定一个bit id,可以推出它对应的特征类型和距离,可以方便地对结果进行解释。

1

2.参数设置

生成药效团指纹首先要计算分子中所有的药效团(化学特征),化学特征的计算需要通过特征工厂实现。

  • 构建特征工厂:
>>> featFactory = ChemicalFeatures.BuildFeatureFactory(fdefName)

使用特征工厂再来构建指纹工厂signature factory用于设置指纹参数。

  • 构建指纹工厂:SigFactory(featFactory, useCounts=False, minPointCount=2, maxPointCount=3, ...)
    featFactory:特征工厂。
    useCounts:默认False。False不考虑指纹频数,并生成SparseBitVect。
    minPointCount:默认为2。生成指纹时包括的最少的药效团数量。
    maxPointCount:默认为3。生成指纹时包括的最多的药效团数量。
>>> from rdkit.Chem.Pharm2D.SigFactory import SigFactory
>>> sigFactory = SigFactory(featFactory, minPointCount=2, maxPointCount=3)
  • 对拓扑距离进行分段:SetBins()
  • 每次修改参数后,都要初始化一下:Init()
  • 计算一下的指纹长度:GetSigSize()
>>> sigFactory.SetBins([(0,2), (2,5), (5,8)])
>>> sigFactory.Init()
>>> sigFactory.GetSigSize()
2988

3.生成2D药效团指纹

指纹工厂中的参数设置完毕,接下来可以生成2D指纹了

  • 计算2D药效团指纹:Gen2DFingerprint(mol, sigFactory, bitInfo, ...)
    mol:要计算指纹的mol对象
    sigFactory:设置了参数的指纹工厂
    bitInfo:获取指纹id及对应原子
>>> from rdkit.Chem.Pharm2D import Generate
>>> mol = Chem.MolFromSmiles('OCC(=O)CCCN')
>>> fp = Generate.Gen2DFingerprint(mol,sigFactory)
>>> print(len(fp))
>>> print(fp.GetNumOnBits())
2988
23
  • 关于指纹每一位所代表特征的信息、特征的距离矩阵等信息,都可以通过signature factory来查看
>>> list(fp.GetOnBits())[:5]
[1, 7, 8, 10, 20]
>>> sigFactory.GetBitDescription(1)
'Acceptor Acceptor |0 1|1 0|'

4.修改FDef设置

如果不想要某个特征,可以直接通过signature factory来跳过某个特征,而不用去修改FDef文件。

  • 查看现有的药效团(特征家族)列表:GetFeatureFamilies()
>>> featFactory.GetFeatureFamilies()
('Donor',
 'Acceptor',
 'NegIonizable',
 'PosIonizable',
 'ZnBinder',
 'Aromatic',
 'Hydrophobe',
 'LumpedHydrophobe')
  • 跳过某个药效团(特征家族):sigFactory.skipFeats
  • 每次修改都要初始化一下:Inint()
  • 再查看一下指纹长度:GetSigSize()
>>> sigFactory.skipFeats=['PosIonizable']
>>> sigFactory.Init()
>>> sigFactory.GetSigSize()
2100
  • 重新生成新的指纹:Gen2DFingerprint()
>>> fp2 = Generate.Gen2DFingerprint(mol,sigFactory)
>>> fp2.GetNumOnBits()
15

5.Gobbi 2D药效团指纹

Rdkit中还有一种用于生成2D药效团指纹的特征定义方式,根据Gobbi等人的设计实现,在 rdkit.Chem.Pharm2D.Gobbi_Pharm2D下有一个预定义的signature factory,已经包含了这些指纹类型,可以直接调用,操作方法类似。

>>> from rdkit.Chem.Pharm2D import Gobbi_Pharm2D, Generate
>>> m = Chem.MolFromSmiles('OCC=CC(=O)O')
>>> fp = Generate.Gen2DFingerprint(m, Gobbi_Pharm2D.factory)
>>> fp.GetNumOnBits()
8
>>> list(fp.GetOnBits())
[23, 30, 150, 154, 157, 185, 28878, 30184]
>>> Gobbi_Pharm2D.factory.GetBitDescription(157)
'HA HD |0 3|3 0|'

本文参考自rdkit官方文档
代码及源文件在这里

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 218,036评论 6 506
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 93,046评论 3 395
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 164,411评论 0 354
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 58,622评论 1 293
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 67,661评论 6 392
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 51,521评论 1 304
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 40,288评论 3 418
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 39,200评论 0 276
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 45,644评论 1 314
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,837评论 3 336
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,953评论 1 348
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,673评论 5 346
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 41,281评论 3 329
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,889评论 0 22
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 33,011评论 1 269
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 48,119评论 3 370
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,901评论 2 355

推荐阅读更多精彩内容