Performance of machine-learning scoring functions in structure-based virtual screening.
Wójcikowski, M., Ballester, P.J., and Siedlecki, P., Sci.Rep. 2017(7)46710
pubmed
こちらの文献は前回使用したODDT(Open Drug Dicovery Toolkit)と同一の筆頭著者によるもので、どうやらドッキングのスコアを機械学習を用いて改善したという内容のようです。しかもODDT同様に、このスコア関数（RF-ScoreVS）を利用できるとのこと！！
文献中では機械学習スコア関数の例としてRF-Score、NNScore、SFCscoreといった名前が挙げられています。このうち、RF-Scoreは今回の文献の著者の一人、De. Pedro Ballester（CRCM）が筆頭著者として2010年に報告しているものです*4。RFはRandom Forestの略のようですが、こちらをversion 1 （RF-Score v1）として、RF-Score v2*5が報告されています。特にRF-Score v3の文献はタイトルが「Improving AutoDock Vina Using Random Forest: The Growing Accuracy of Binding Affinity Prediction by the Effective Exploitation of Larger Data Sets.」、とAutoDock Vinaを使った以上とても気になる文献です（アクセスできなかった・・・）。
x = Rescored_df['RFScoreVS_v2']
y = Rescored_df['vina_affinity']
plt.scatter(x,y)
plt.title("RF-Score-VS vs AutoDock Vina")
plt.xlabel("RF-Score-VS")
plt.ylabel("AutoDock Vina")
plt.grid(True)
*3:Machine-learning scoring functions to improve structure-based binding affinity prediction and virtual screening
Ain, Q.U., Aleksandrova, A., Roessler, F.D., and Ballester P.J. WIREs Comput. Mol. Sci.2015(5)405
*4:A machine learning approach to predicting protein-ligand binding affinity with applications to molecular docking. Ballester, P.J. and Mitchell., J.B.O., Bioinformatics. 2010(26)1169
*5:Does a more precise chemical description of protein-ligand complexes lead to more accurate prediction of binding affinity? Ballester, P.J., Schreyer, A., Blundell. T.L., J. Chem. Inf. Model. 2014(54)944)、RF-Score v3 ((Li, H., Leung, KS., Wong, M.H., Ballester P.J., Mol. Inform. 2015(34)115
from rdkit.Chem.FeatMaps import FeatMaps
from rdkit.Chem import rdMolAlign
defFeatMapScorer(mol, tmp):
#比較対象(tmp)との間でアラインメントをとる
o3d = rdMolAlign.GetO3A(mol,tmp)
o3d.Align()
#feature mapのパラメーターの辞書を作成（デフォルトの値を使用）
fmParams = {}
for k in Ffactory.GetFeatureFamilies():
fparams = FeatMaps.FeatMapParams()
fmParams[k] = fparams
#featureのうちより一般的なものだけを残すための設定
keep = ('Donor','Acceptor','NegIonizable','PosIonizable','Aromatic')
featLists = []
for m in [mol, tmp]:
#featureをまず全部取得
rawFeats = fdef.GetFeaturesForMol(m)
#keepにある欲しいものだけリストに追加
featLists.append([f for f in rawFeats if f.GetFamily() in keep])
#FeatMap オブジェクトを作成
fms = [FeatMaps.FeatMap(feats=x, weights=[1]*len(x), params=fmParams) for x in featLists]
#feature mapに対するスコアを取得
score = fms[0].ScoreFeats(featLists[1]) / min(fms[0].GetNumFeatures(), len(featLists[1]))
#scoreを返すreturn score
#ODDTでSDF書き出し
DockAllResultOutput = oddt.toolkit.Outputfile('sdf', 'ODDT_Dock_All_Result.sdf')
for mol in All_Result:
DockAllResultOutput.write(mol)
DockAllResultOutput.close()
#RDKitのPandasToolsで読み込み
All_Dock_df = PandasTools.LoadSDF('./ODDT_Dock_All_Result.sdf')
"""RDKit ERROR: [22:07:38] Explicit valence for atom # 2 C, 5, is greater than permittedRDKit ERROR: [22:07:38] ERROR: Could not sanitize molecule ending on line 42625RDKit ERROR: [22:07:38] ERROR: Explicit valence for atom # 2 C, 5, is greater than permitted・・・以下エラーが続く"""
*11:Wójcikowski, M., Zielenkiewicz, P., and Siedlecki, P. Open Drug Discovery Toolkit (ODDT): a new open-source player in the drug discovery field Journal of Cheminformatics 2015(7)26
from rdkit import Chem
from rdkit.Chem import AllChem
ortho_biphenyl= Chem.MolFromSmarts('c1ccc(c(*)c1)c1ccccc1')
P_Matched2_suppl = Chem.ForwardSDMolSupplier('./P_Matched2.sdf')
P_Matched2_mols = [mol for mol in P_Matched2_suppl if mol isnotNone]
print(len(P_Matched2_mols))
#997#オルト置換ありのリスト
ortho_cpds = []
#オルト置換なしのリスト
w_o_ortho = []
for mol in P_Matched2_mols:
if mol.HasSubstructMatch(ortho_biphenyl):
ortho_cpds.append(mol)
else:
w_o_ortho.append(mol)
print('Number of ortho substituted compounds', len(ortho_cpds))
# Number of ortho substituted compounds 494print('Number of ortho unsubstituted compounds', len(w_o_ortho))
# Number of ortho unsubstituted compounds 503
maccs_fps = [AllChem.GetMACCSKeysFingerprint(mol) for mol in fp_mols]
rdkit_fps = [AllChem.RDKFingerprint(mol) for mol in fp_mols]
ECFP_fps = [AllChem.GetMorganFingerprintAsBitVect(mol, 2, 2048) for mol in fp_mols]
FCFP_fps = [AllChem.GetMorganFingerprintAsBitVect(mol, 2, 2048, useFeatures=True) for mol in fp_mols]
Apair_fps = [AllChem.GetAtomPairFingerprint(mol) for mol in fp_mols]
from rdkit.Avalon import pyAvalonTools
avalon_fps = [pyAvalonTools.GetAvalonFP(mol) for mol in fp_mols]
defcommon_bit(fps):
fp_size = len(fp2arr(fps[0]))
sumarr = np.zeros(fp_size)
for i in fps:
arr = fp2arr(i)
sumarr += arr
mean = sumarr.mean()
threshold = math.floor(mean)
common_bit = []
for x in sumarr:
if x >= threshold:
common_bit.append(1)
else:
common_bit.append(0)
return fp_size, threshold, sumarr, common_bit
import statistics
Scores_STD = []
for l in Scores:
l_mean = statistics.mean(l)
l_stdev = statistics.stdev(l)
s = [(i-l_mean) / l_stdev for i in l]
Scores_STD.append(s)
標準化したスコアを先と同様にプロットします。
for i inrange(6):
plt.subplot(2, 3, i+1)
plt.hist(Scores_STD[i], bins=20)
plt.title(fp_names[i])
Integrated_SC_arr = np.zeros(494)
for i in Scores_STD:
Integrated_SC_arr += np.array(i)
#最小値を求める
int_min = min(Integrated_SC_arr)
#最小値の絶対値分シフト
Integrated_SC = [i + abs(int_min) for i in Integrated_SC_arr]
#グラフ化
plt.hist(Integrated_SC, bins=20)
当てずっぽうはさておき、LipinskiのRule of 5からは真っ先に弾かれてしまいそうなこれらの分子が経口薬として上市に至っている、という成功例がPD-1/PD-L1 結合阻害剤においても巨大化、対称化の流れを進める背景にもあったりするのかな、と想像してしまいます。創薬の現場からみたらどうなんでしょう？？？
from rdkit import Chem
aligned_lig_sup = Chem.SDMolSupplier('./aligned_ligands.sdf')
aligned_ligs = [x for x in aligned_lig_sup if x isnotNone]
以下のようなエラーがたくさん出てきました。
RDKit ERROR: [20:29:00] ERROR: non-ring atom 14 marked aromatic
RDKit ERROR: [21:10:06] non-ring atom 4 marked aromatic
RDKit ERROR: [21:10:06] ERROR: Could not sanitize molecule ending on line 72
for x in nbrs:
print('AtomID: ', x.GetIdx(), 'Symbol: ', x.GetSymbol(), 'CIPRank: ', x.GetProp('_CIPRank'))
#AtomID: 30 Symbol: C CIPRank: 4#AtomID: 32 Symbol: C CIPRank: 39#AtomID: 26 Symbol: N CIPRank: 42
chiral_centerH = test_m3H.GetAtomWithIdx(31)
nbrsH = chiral_centerH.GetNeighbors()
for i in nbrsH:
print('AtomID: ', i.GetIdx(), 'Symbol: ', i.GetSymbol())
#AtomID: 30 Symbol: C
#AtomID: 32 Symbol: C
#AtomID: 26 Symbol: N
#AtomID: 79 Symbol: H
defPSearch3(mol):
#水素原子を付加しない
smi = Chem.MolToSmiles(mol, isomericSmiles=False)
mol_re = Chem.MolFromSmiles(smi)
AllChem.Compute2DCoords(mol_re)
#フィーチャーをそもそも持っているか？
match, mList = EmbedLib.MatchPharmacophoreToMol(mol_re, Ffactory, pcophore)
if match == True:
#距離の検証#距離行列の取得
bounds = rdDistGeom.GetMoleculeBoundsMatrix(mol_re)
#ファーマコフォアのマッチング
pList = EmbedLib.GetAllPharmacophoreMatches(mList,bounds,pcophore)
#pListのそれぞれについてFeatureとマッチする原子のidを取得する
num_match = len(pList)
phMatches = []
for i inrange(num_match):
num_feature = len(pList[i])
phMatch = []
for j inrange(num_feature):
phMatch.append(pList[i][j].GetAtomIds())
phMatches.append(phMatch)
#原子のidとファーマコフォアをもとに3次元構造を埋め込む
Generated_embeds = []
for phMatch in phMatches:
bm,embeds,nFail =EmbedLib.EmbedPharmacophore(mol_re, phMatch, pcophore,count=5, silent=1)
Generated_embeds.append(len(embeds))
returnsum(Generated_embeds)
else:
return0
関数3の検証
活性化合物群に適用します。
PScore3_list = []
for x in chain_mols:
num_embeds = PSearch3(x)
PScore3_list.append(num_embeds)
#Molオブジェクトのプロパティに結果を保持させる
x.SetIntProp('PScore', num_embeds)
biphenyl_suppl = Chem.SDMolSupplier('./biphenyl_library.sdf')
biphenyl_mols = [x for x in biphenyl_suppl if x isnotNone]
biphenyl_PScores = []
for x in biphenyl_mols:
num_embeds = PSearch3(x)
biphenyl_PScores.append(num_embeds)
x.SetIntProp('PScore', num_embeds)
ヒット（一つ以上ファーマコフォアを満たす構造がembedできた化合物）の数を確認します。
Hit_list = [i for i, x inenumerate(biphenyl_PScores) if x != 0]
print(len(Hit_list))
# 2463
合計2463個の化合物が残りました。約半数に絞り込めたようです。
上記の計算過程で、
RDKit ERROR: [23:37:56] DEBUG: Removed embedding due to c[00:09:28] UFFTYPER: Unrecognized charge state for atom: 1
print([i for i, x inenumerate(biphenyl_PScores) if x == max(biphenyl_PScores)])
# [1351]
max_mol = biphenyl_mols[1351]
max_mol.GetProp('PScore')
#'99'
Draw.MolToImage(max_mol)
数は99で、構造は以下でした。
・・・ほとんどsp2です。溶解度悪そう、、、
とりあえず結果をSDFで出力しておきます。
#保持するためのpropertyを確認
check_mol = biphenyl_mols[0]
prop_names = [p for p in check_mol.GetPropNames()]
print(prop_names)
# ['original_id', 'MW', 'MolLogP', 'NumHDonors', 'NumHAcceptors', 'NumRotatableBonds', 'TPSA', 'PScore']
writer = Chem.SDWriter('PSearched.sdf')
writer.SetProps(prop_names)
for mol in biphenyl_mols:
writer.write(mol)
writer.close()
import os
from rdkit import RDConfig
from rdkit.Chem import AllChem
from rdkit import Chem
chain_suppl = Chem.SDMolSupplier('./chain_compounds.sdf', removeHs=False)
chain_mols = [x for x in chain_suppl if x isnotNone]
# コピペfrom rdkit.Chem.Features import ShowFeats
defshowMolFeats(mol,factory,viewer,name="mol",showOnly=True):
featLabel = f'{name}-feats'
dirLabel = f'{name}-feats-dirs'if(showOnly):
ShowFeats._globalSphereCGO = []
ShowFeats._globalArrowCGO = []
# workaround for what is either a bug in the way pymol handles CGOs# or a gap in my understanding of how it should work
viewer.server.resetCGO(featLabel)
viewer.server.resetCGO(dirLabel)
viewer.server.sphere((0, 0, 0), .01, (1, 0, 1), featLabel)
viewer.server.cylinder((0, 0, 0), (.01, .01, .01), .01, (1, 0, 1), dirLabel)
ShowFeats.ShowMolFeats(mol,factory,v,name=name,featLabel=featLabel,dirLabel=dirLabel,
useFeatDirs=False,showOnly=showOnly)
viewer.server.renderCGO(ShowFeats._globalSphereCGO,featLabel,1)
viewer.server.renderCGO(ShowFeats._globalArrowCGO,dirLabel,1)
suppl_8YQ = Chem.SDMolSupplier('./8YQ_noH.sdf')
mol_8YQ = [x for x in suppl_8YQ if x isnotNone][0]
showMolFeats(mol_8YQ,fdef,v)
png=v.GetPNG()
display.display(png)
