分子基礎(chǔ)操作與藥效團(tuán)查找
文章目錄
- 1.原子操作
- 2.鍵操作
- 3.環(huán)操作
- 4.手動(dòng)實(shí)現(xiàn)氧族藥效團(tuán)查找
1.原子操作
在rdkit中，分子中的每一個(gè)原子都是對(duì)象措伐，可以通過(guò)原子對(duì)象的屬性和函數(shù)來(lái)獲取各種信息。
- 對(duì)原子進(jìn)行遍歷：m.GetAtoms()
- 獲取原子索引：GetIdx()
- 獲取原子序號(hào)：GetAtomicNum()
- 獲取原子符號(hào)：GetSymbol()
- 獲取原子連接數(shù)（受H是否隱藏影響）：GetDegree()
- 獲取原子總連接數(shù)（與H是否隱藏?zé)o關(guān)）：GetTotalDegree()
- 獲取原子形式電荷：GetFormalCharge()
- 獲取原子雜化方式：GetHybridization()
- 獲取原子顯式化合價(jià)：GetExplicitValence()
- 獲取原子隱式化合價(jià)：GetImplicitValence()
- 獲取原子總的化合價(jià)：GetTotalValence()
- ...
>>> from rdkit import Chem
>>> m = Chem.MolFromSmiles('C1OC1')
>>> print('\t'.join(['id', 'num', 'symbol', 'degree', 'charge', 'hybrid']))
>>> for atom in m.GetAtoms():
>>> print(atom.GetIdx(), end='\t')
>>> print(atom.GetAtomicNum(), end='\t')
>>> print(atom.GetSymbol(), end='\t')
>>> print(atom.GetDegree(), end='\t')
>>> print(atom.GetFormalCharge(), end='\t')
>>> print(atom.GetHybridization())
id num symbol degree charge hybrid
0 6 C 2 0 SP3
1 8 O 2 0 SP3
2 6 C 2 0 SP3
- 也可以通過(guò)索引獲取原子：GetAtomWithIdx()
>>> print(m.GetAtomWithIdx(0).GetSymbol())
C
- 獲取相連的原子：GetNeighbors()
>>> atom = m.GetAtomWithIdx(1)
>>> [x.GetAtomicNum() for x in atom.GetNeighbors()]
[6, 6]
2.鍵操作
同樣尾序，每一個(gè)鍵也都是對(duì)象挠轴，可以通過(guò)屬性和函數(shù)來(lái)獲取鍵的信息。
- 對(duì)鍵進(jìn)行遍歷：m.GetBonds()
- 獲取鍵的索引：GetIdx()
- 獲取鍵的類(lèi)型：GetGetBondType()
- 以數(shù)字形式顯示鍵的類(lèi)型：GetBondTypeAsDouble()
- 是否為芳香鍵：GetIsAromatic()
- 是否為共軛鍵：GetIsConjugated()
- 是否在環(huán)中：IsInRing()
- 是否在n元環(huán)中：IsInRingSize(n)
- 獲取起始原子：GetBeginAtom()
- 獲取末尾原子：GetEndAtom()
- ...
>>> m = Chem.MolFromSmiles('C=C-C=C')
>>> print('\t'.join(['id', 'type', 'double', 'aromic', 'conjug', 'ring', 'begin', 'end']))
>>> for bond in m.GetBonds():
>>> print(bond.GetIdx(), end='\t')
>>> print(bond.GetBondType(), end='\t')
>>> print(bond.GetBondTypeAsDouble(), end='\t')
>>> print(bond.GetIsAromatic(), end='\t')
>>> print(bond.GetIsConjugated(), end='\t')
>>> print(bond.IsInRing(), end='\t')
>>> print(bond.GetBeginAtomIdx(), end='\t')
>>> print(bond.GetEndAtomIdx())
id type double aromic conjug ring begin end
0 DOUBLE 2.0 False True False 0 1
1 SINGLE 1.0 False True False 1 2
2 DOUBLE 2.0 False True False 2 3
- 也可以通過(guò)索引獲取鍵：GetBondWithIdx()
>>> print(m.GetBondWithIdx(0).GetBondType())
DOUBLE
3.環(huán)操作
- 查看原子的環(huán)信息
>>> m = Chem.MolFromSmiles('OC1C2C1CC2')
>>> print(m.GetAtomWithIdx(0).IsInRing())
False
>>> print(m.GetAtomWithIdx(2).IsInRing())
True
>>> print(m.GetAtomWithIdx(2).IsInRingSize(3))
True
>>> m
atom_bond_ring_1.png
- 但是要注意翩肌，IsInRingSize()函數(shù)只能判斷是否在最小的環(huán)中
>>> print(m.GetAtomWithIdx(2).IsInRingSize(3))
True
>>> print(m.GetAtomWithIdx(2).IsInRingSize(5))
False
- 可以查看所有最小環(huán)(smallest set of smallest rings, SSSR）的信息：GetSymmSSSR()
>>> ssr = Chem.GetSymmSSSR(m)
>>> print(len(ssr))
>>> print(list(ssr[0]))
>>> print(list(ssr[1]))
2
[1, 2, 3]
[4, 5, 2, 3]
- 直接獲取環(huán)的信息：GetRingInfo()
- 查看一共有幾個(gè)環(huán)：NumRings()
- 查看原子在幾個(gè)環(huán)中：NumAtomRings()
- 查看id為n的原子是否在n1元環(huán)中.IsAtomInRingOfSize(n, n1)
- 查看id為n的鍵是否在n1元環(huán)中.IsAtomInRingOfSize(n , n1)
>>> ri = m.GetRingInfo()
>>> print(ri.NumRings())
2
>>> print(ri.NumAtomRings(6))
0
>>> print(ri.IsAtomInRingOfSize(1,3))
True
>>> print(ri.IsBondInRingOfSize(1,3))
True
4.手動(dòng)實(shí)現(xiàn)氧族藥效團(tuán)查找
這里只做一個(gè)簡(jiǎn)單的演示模暗，通過(guò)原子屬性查找目標(biāo)藥效團(tuán)，更復(fù)雜的操作類(lèi)似念祭。
- 假設(shè)要查找的氧族氫供體標(biāo)準(zhǔn)：氧或硫原子兑宇，不帶電荷，含有1個(gè)氫粱坤。
- 氧族氫受體標(biāo)準(zhǔn)（部分）：不帶氫原子隶糕，化合價(jià)為2瓷产，且不與氮原子相連。
>>> def pharmacophore(m):
>>> atomPharma = {}
>>> # 定義氧族原子
>>> chalcogen = [8, 16]
>>> mol = Chem.MolFromSmiles(m)
>>> mol = Chem.AddHs(mol)
>>> # 開(kāi)始查找
>>> for atom in mol.GetAtoms():
>>> p = []
>>> if atom.GetAtomicNum() == 1 or atom.GetAtomicNum() not in chalcogen:
>>> continue
>>> else:
>>> # 氫供體
>>> if atom.GetFormalCharge() == 0:
>>> nbrs = [x for x in atom.GetNeighbors()]
>>> HDflag = False
>>> for nbr in nbrs:
>>> if nbr.GetAtomicNum() == 1:
>>> HDflag = True
>>> if HDflag:
>>> p.append('HDonor')
>>> # 氫受體
>>> if atom.GetTotalValence() == 2:
>>> nbrs = [x for x in atom.GetNeighbors()]
>>> HAflag_1 = True
>>> HAflag_2 = True
>>> if len(nbrs) == 1:
>>> nbr = nbrs[0]
>>> if nbr.GetAtomicNum() == 7:
>>> HAflag_1 = False
>>> else:
>>> for nbr in nbrs:
>>> if nbr.GetAtomicNum() == 1:
>>> HAflag_2 = False
>>> if HAflag_1 and HAflag_2:
>>> p.append('HAcceptor')
>>> atomPharma[atom.GetIdx()] = [atom.GetAtomicNum(), ' '.join(p)]
>>> return atomPharma
>>>
>>> m = 'COC(=O)O'
>>> res = pharmacophore(m)
>>> res
{1: [8, 'HAcceptor'], 3: [8, 'HAcceptor'], 4: [8, 'HDonor']}
更簡(jiǎn)單的操作枚驻，直接寫(xiě)成SMARTS查找就可以了
- 氫供體：[O,S;H1;+0]
- 氫受體（部分）：[O;H0;v2;!$(O=N-*)]
>>> def pharmacophore_smarts(m):
>>> # 定義smarts
>>> HDonorSmarts = '[O,S;H1;+0]'
>>> HAcceptorSmarts = '[O;H0;v2;!$(O=N-*)]'
>>> HDonor = Chem.MolFromSmarts(HDonorSmarts)
>>> HAcceptor = Chem.MolFromSmarts(HAcceptorSmarts)
>>> atomPharma = {}
>>> mol = Chem.MolFromSmiles(m)
>>> # 氫供體
>>> HDonors = mol.GetSubstructMatches(HDonor)
>>> for i in HDonors:
>>> atom = mol.GetAtomWithIdx(i[0])
>>> atomPharma[atom.GetIdx()] = [atom.GetAtomicNum(), 'HDonor']
>>> # 氫受體
>>> HAcceptors = mol.GetSubstructMatches(HAcceptor)
>>> for i in HAcceptors:
>>> atom = mol.GetAtomWithIdx(i[0])
>>> atom_prop = atomPharma.get(atom.GetIdx(), [])
>>> if atom_prop:
>>> atom_prop[1] = atom_prop[1] + ' HAcceptor'
>>> atomPharma.update({atom.GetIdx():atom_prop})
>>> else:
>>> atomPharma[atom.GetIdx()] = [atom.GetAtomicNum(), 'HAcceptor']
>>> return atomPharma
>>>
>>> res = pharmacophore_smarts(m)
>>> res
{4: [8, 'HDonor'], 1: [8, 'HAcceptor'], 3: [8, 'HAcceptor']}
