• 一、分子片段生成
• 二、片段指紋生成
• 三、指紋重要性分析

一、分子片段生成

分子片段（Molecular Fragments）是一組相連的原子，并可能包含有相關(guān)官能團(tuán)。在rdkit中提供了一系列用于分析、操作分子片段的工具。說起來比較抽象，操作起來也比較抽象。

• 獲取官能團(tuán)庫：RDConfig.RDDataDir目錄下的'FunctionalGroups.txt'
• 根據(jù)官能團(tuán)庫實(shí)例化一個(gè)參數(shù)器：FragmentCatalog.FragCatParams()

>>> import os
>>> from rdkit import Chem
>>> from rdkit.Chem import Draw
>>> from rdkit.Chem import RDConfig
>>> from rdkit.Chem import FragmentCatalog
>>> fName = os.path.join(RDConfig.RDDataDir, 'FunctionalGroups.txt')
>>> fparams = FragmentCatalog.FragCatParams(1, 6, fName)

• 查看庫中包含的官能團(tuán)數(shù)量：GetNumFuncGroups()
• 查看每個(gè)官能團(tuán)對(duì)應(yīng)的基團(tuán)：GetFuncGroup()

>>> print(fparams.GetNumFuncGroups())
>>> m = fparams.GetFuncGroup(0)
>>> m
39

1

• 傳入?yún)?shù)器，創(chuàng)建一個(gè)片段存儲(chǔ)器：FragmentCatalog.FragCatalog()
  產(chǎn)生的分子片段都會(huì)存儲(chǔ)在該對(duì)象中
• 創(chuàng)建一個(gè)片段生成器：FragmentCatalog.FragCatGenerator()
  通過該對(duì)象生成片段
• 計(jì)算分子片段：AddFragsFromMol()
• 查看分子片段數(shù)量：GetNumEntries()

>>> m = Chem.MolFromSmiles('OCC=CC(=O)O')
>>> fcat = FragmentCatalog.FragCatalog(fparams)
>>> fcgen = FragmentCatalog.FragCatGenerator()
>>> fcgen.AddFragsFromMol(m, fcat)
>>> fcat.GetNumEntries()
3

• 通過存儲(chǔ)器查看片段：fcat.GetEntryDescription()
  尖括號(hào)中的內(nèi)容：表示與片段相連的官能團(tuán)，以下面的結(jié)果為例：
  第0號(hào)片段中，對(duì)應(yīng)著一個(gè)乙基片段，該乙基與一個(gè)羥基相連。
  第1號(hào)片段中，對(duì)應(yīng)著一個(gè)乙烯片段，該乙烯與一個(gè)羧基相連。

>>> print(fcat.GetEntryDescription(0))
>>> print(fcat.GetEntryDescription(1))
C<-O>C
C=C<-C(=O)O>

關(guān)于官能團(tuán)的詳細(xì)信息，可以通過下述方法獲?。?/p>

• 向存儲(chǔ)器傳入分子片段id，獲取片段中所包含的官能團(tuán)編號(hào)：fcat.GetEntryFuncGroupIds()
• 向參數(shù)器傳入官能團(tuán)編號(hào)，獲取官能團(tuán)對(duì)應(yīng)的mol對(duì)象：fparams.GetFuncGroup()

>>> print(list(fcat.GetEntryFuncGroupIds(0)))
>>> funcgroup = fparams.GetFuncGroup(34)
>>> print(Chem.MolToSmarts(funcgroup))
>>> print(funcgroup.GetProp('_Name'))
[34]
*-[O&D1]
-O

提取得到的片段是層級(jí)結(jié)構(gòu)，小片段在最底層，逐漸合并形成大片段?？梢圆榭匆粋€(gè)小片段形成了哪些大片段。

• 根據(jù)id獲取片段：fcat.GetEntryDescription()
• 獲取上級(jí)片段id：fcat.GetEntryDownIds()

>>> print(fcat.GetEntryDescription(0))
>>> list(fcat.GetEntryDownIds(0))
C<-O>C
[2]

• 根據(jù)上級(jí)片段id，獲取上級(jí)片段信息

>>> fcat.GetEntryDescription(2)
'C<-C(=O)O>=CC<-O>'

二、片段指紋生成

• 先將多個(gè)分子的片段匯總到一個(gè)片段存儲(chǔ)器中

>>> ms = [Chem.MolFromSmiles('OCC(NC1CC1)CCC'), Chem.MolFromSmiles('OCC=CC(=O)O')]
>>> fcat = FragmentCatalog.FragCatalog(fparams)
>>> for m in ms:
>>>     fcgen.AddFragsFromMol(m, fcat)
>>> fcat.GetNumEntries()
17

存儲(chǔ)器收集完所有片段后，再用它來生成分子指紋

• 創(chuàng)建一個(gè)片段指紋生成器：FragFPGenerator()
• 傳入分子和存儲(chǔ)器來生成指紋：GetFPForMol(mol, fcat)
• 以字符串形式查看指紋：ToBitString()
• 查看指紋中哪些位是有效的：GetOnBits()

>>> fpgen = FragmentCatalog.FragFPGenerator()
>>> fp1 = fpgen.GetFPForMol(ms[1], fcat)
>>> print(fp1.ToBitString())
>>> print(list(fp1.GetOnBits()))
10000000000000011
[0, 15, 16]

可以用處理一般分子指紋的方法來處理片段分子指紋，例如尋找相同的片段

• 先對(duì)分子指紋做一步“&”位運(yùn)算，兩個(gè)指紋都為1時(shí)，結(jié)果為1，否則為0
• 獲取兩個(gè)指紋中都出現(xiàn)的片段：GetOnBits()
• 查看片段信息：GetEntryDescription()

>>> fp0 = fpgen.GetFPForMol(ms[0], fcat)
>>> andfp = fp0 & fp1
>>> onbit = list(andfp.GetOnBits())
>>> fcat.GetEntryDescription(onbit[0])
'C<-O>C'

也可以按上述思路查看一下哪些片段導(dǎo)致了分子的不同

• 對(duì)分子指紋做“^”運(yùn)算，兩個(gè)指紋相同時(shí)，結(jié)果為0，否則為1。再做一步“&”運(yùn)算
• 按上述方法查看相異片段

>>> dis = fp0 ^ fp1
>>> combinedfp = fp0 & dis
>>> onbit = list(combinedfp.GetOnBits())
>>> fcat.GetEntryDescription(onbit[-1])
'CCCC(C<-O>)N<-cPropyl>'

三、指紋重要性分析

這里主要介紹指紋對(duì)離散標(biāo)簽的重要性分析。在rdkit.ML.InfoTheory.rdInfoTheory.InfoBitRanker中提供了對(duì)指紋分析的功能。這個(gè)類可以根據(jù)分子指紋和離散標(biāo)簽，對(duì)特征進(jìn)行計(jì)算和排序，看看哪些特征對(duì)活性比較重要。

• 先對(duì)163個(gè)分子生成片段指紋，完整代碼如下

>>> suppl = Chem.SDMolSupplier('data/bzr.sdf')
>>> sdms = [x for x in suppl]

>>> # 獲取官能團(tuán)庫
>>> fName = os.path.join(RDConfig.RDDataDir, 'FunctionalGroups.txt')
>>> # 片段參數(shù)器
>>> fparams = FragmentCatalog.FragCatParams(1, 6, fName)
>>> # 片段存儲(chǔ)器
>>> fcat = FragmentCatalog.FragCatalog(fparams)
>>> # 片段生成器
>>> fcgen = FragmentCatalog.FragCatGenerator()
>>> # 片段指紋生成器
>>> fpgen = FragmentCatalog.FragFPGenerator()
>>> # 匯總所有片段
>>> for m in sdms:
>>>     fcgen.AddFragsFromMol(m, fcat)
>>> # 生成片段指紋
>>> fps = [fpgen.GetFPForMol(x, fcat) for x in sdms]
>>> print(len(fps), fps[0].GetNumBits())
163 8266

• 信息增益（infoGain）分析，實(shí)例化一個(gè)排序?qū)ο螅篒nfoBitRanker(nBits, nClasses, infoType)
  nBits：指紋長(zhǎng)度
  nClasses：類別數(shù)量，需要和標(biāo)簽滿足的關(guān)系：0 <= 標(biāo)簽 < 類別數(shù)量
  infoType：度量指標(biāo)。默認(rèn)使用rdInfoTheory.InfoType.ENTROPY，即信息增益作為比較標(biāo)準(zhǔn)，它反映了使用某個(gè)特征進(jìn)行分類后，系統(tǒng)混亂程度降低的多少，數(shù)值越大表明特征越重要。

>>> from rdkit.ML import InfoTheory
>>> ranker = InfoTheory.InfoBitRanker(len(fps[0]), 2)

• 獲取每個(gè)分子的活性信息：GetDoubleProp('ACTIVITY')
• 以7作為標(biāo)準(zhǔn)對(duì)活性離散化，大于7為1，小于7為0
• 根據(jù)指紋和類別進(jìn)行投票：AccumulateVotes(fp, act)
• 獲取前5個(gè)重要特征：GetTopN(5)
• 依次輸出特征id、信息增益、特征為1類別中的無活性分子數(shù)、特征為1類別中的有活性分子數(shù)。

>>> acts = [x.GetDoubleProp('ACTIVITY') for x in sdms]
>>> for i,fp in enumerate(fps):
>>>     act = int(acts[i]>7)
>>>     ranker.AccumulateVotes(fp,act)
>>> top5 = ranker.GetTopN(5)
>>> for id, gain, n0, n1 in top5:
>>>     print(int(id), '%.3f'%gain, int(n0), int(n1))
698 0.081 20 17
222 0.073 23 25
378 0.073 30 43
196 0.073 30 43
1207 0.073 0 25

• 加入偏置，以信息增益為例，重新設(shè)置infoType
• 設(shè)置偏置類別：SetBiasList()
  在這種模式下，一個(gè)特征與所設(shè)置了偏置類別的相關(guān)性要高于所有非偏置類別，例如設(shè)置偏置類別為4，某位特征為1對(duì)應(yīng)的標(biāo)簽中，類別為4的數(shù)量應(yīng)該大于其他類別的數(shù)量。

>>> ranker = InfoTheory.InfoBitRanker(len(fps[0]), 2, InfoTheory.InfoType.BIASENTROPY)
>>> ranker.SetBiasList((0,))
>>> acts = [x.GetDoubleProp('ACTIVITY') for x in sdms]
>>> for i,fp in enumerate(fps):
>>>     act = 0 if acts[i]<7 else 1
>>>     ranker.AccumulateVotes(fp, act)
>>> top5 = ranker.GetTopN(5)
>>> for id, gain, n0, n1 in top5:
>>>     print(int(id), '%.3f'%gain, int(n0), int(n1))
698 0.081 20 17
222 0.073 23 25
378 0.073 30 43
196 0.073 30 43
2375 0.062 5 0

• 使用卡方檢驗(yàn)（chi squared test），將infoType設(shè)置為如下參數(shù)，其他相同

>>> ranker = InfoTheory.InfoBitRanker(len(fps[0]), 2, InfoTheory.InfoType.CHISQUARE)
>>> for i,fp in enumerate(fps):
>>>     act = int(acts[i]>7)
>>>     ranker.AccumulateVotes(fp, act)
>>> top5 = ranker.GetTopN(5)
>>> for id, gain, n0, n1 in top5:
>>>     print(int(id), '%.3f'%gain, int(n0), int(n1))
698 20.023 20 17
222 17.451 23 25
378 16.242 30 43
196 16.242 30 43
2375 14.861 5 0

• 帶偏置的卡方檢驗(yàn)同上，重新設(shè)置infoType為InfoTheory.InfoType.BIASCHISQUARE即可

本文參考自rdkit官方文檔。
代碼及源文件在這里。