5
|
1 #!/usr/bin/env python
|
0
|
2 __all__ = ['weights_max', 'weights_mean', 'weights_none', 'default']
|
|
3
|
|
4 # Silicos-it
|
|
5 from errors import WrongArgument
|
|
6
|
|
7 # RDKit
|
|
8 from rdkit.Chem import Descriptors
|
|
9 from rdkit import Chem
|
|
10
|
|
11 # General
|
|
12 from copy import deepcopy
|
|
13 from math import exp, log
|
|
14 import sys, os, re
|
|
15 import argparse
|
|
16
|
|
17 def check_filetype(filepath):
|
|
18 mol = False
|
|
19 for line in open(filepath):
|
|
20 if line.find('$$$$') != -1:
|
|
21 return 'sdf'
|
|
22 elif line.find('@<TRIPOS>MOLECULE') != -1:
|
|
23 return 'mol2'
|
|
24 elif line.find('ligand id') != -1:
|
|
25 return 'drf'
|
|
26 elif re.findall('^InChI=', line):
|
|
27 return 'inchi'
|
|
28 elif re.findall('^M\s+END', line):
|
|
29 mol = True
|
|
30
|
|
31 if mol:
|
|
32 # END can occures before $$$$, so and SDF file will
|
|
33 # be recognised as mol, if you not using this hack'
|
|
34 return 'mol'
|
|
35 return 'smi'
|
|
36
|
|
37 AliphaticRings = Chem.MolFromSmarts('[$([A;R][!a])]')
|
|
38
|
|
39 AcceptorSmarts = [
|
|
40 '[oH0;X2]',
|
|
41 '[OH1;X2;v2]',
|
|
42 '[OH0;X2;v2]',
|
|
43 '[OH0;X1;v2]',
|
|
44 '[O-;X1]',
|
|
45 '[SH0;X2;v2]',
|
|
46 '[SH0;X1;v2]',
|
|
47 '[S-;X1]',
|
|
48 '[nH0;X2]',
|
|
49 '[NH0;X1;v3]',
|
|
50 '[$([N;+0;X3;v3]);!$(N[C,S]=O)]'
|
|
51 ]
|
|
52 Acceptors = []
|
|
53 for hba in AcceptorSmarts:
|
|
54 Acceptors.append(Chem.MolFromSmarts(hba))
|
|
55
|
|
56 StructuralAlertSmarts = [
|
|
57 '*1[O,S,N]*1',
|
|
58 '[S,C](=[O,S])[F,Br,Cl,I]',
|
|
59 '[CX4][Cl,Br,I]',
|
|
60 '[C,c]S(=O)(=O)O[C,c]',
|
|
61 '[$([CH]),$(CC)]#CC(=O)[C,c]',
|
|
62 '[$([CH]),$(CC)]#CC(=O)O[C,c]',
|
|
63 'n[OH]',
|
|
64 '[$([CH]),$(CC)]#CS(=O)(=O)[C,c]',
|
|
65 'C=C(C=O)C=O',
|
|
66 'n1c([F,Cl,Br,I])cccc1',
|
|
67 '[CH1](=O)',
|
|
68 '[O,o][O,o]',
|
|
69 '[C;!R]=[N;!R]',
|
|
70 '[N!R]=[N!R]',
|
|
71 '[#6](=O)[#6](=O)',
|
|
72 '[S,s][S,s]',
|
|
73 '[N,n][NH2]',
|
|
74 'C(=O)N[NH2]',
|
|
75 '[C,c]=S',
|
|
76 '[$([CH2]),$([CH][CX4]),$(C([CX4])[CX4])]=[$([CH2]),$([CH][CX4]),$(C([CX4])[CX4])]',
|
|
77 'C1(=[O,N])C=CC(=[O,N])C=C1',
|
|
78 'C1(=[O,N])C(=[O,N])C=CC=C1',
|
|
79 'a21aa3a(aa1aaaa2)aaaa3',
|
|
80 'a31a(a2a(aa1)aaaa2)aaaa3',
|
|
81 'a1aa2a3a(a1)A=AA=A3=AA=A2',
|
|
82 'c1cc([NH2])ccc1',
|
|
83 '[Hg,Fe,As,Sb,Zn,Se,se,Te,B,Si,Na,Ca,Ge,Ag,Mg,K,Ba,Sr,Be,Ti,Mo,Mn,Ru,Pd,Ni,Cu,Au,Cd,Al,Ga,Sn,Rh,Tl,Bi,Nb,Li,Pb,Hf,Ho]',
|
|
84 'I',
|
|
85 'OS(=O)(=O)[O-]',
|
|
86 '[N+](=O)[O-]',
|
|
87 'C(=O)N[OH]',
|
|
88 'C1NC(=O)NC(=O)1',
|
|
89 '[SH]',
|
|
90 '[S-]',
|
|
91 'c1ccc([Cl,Br,I,F])c([Cl,Br,I,F])c1[Cl,Br,I,F]',
|
|
92 'c1cc([Cl,Br,I,F])cc([Cl,Br,I,F])c1[Cl,Br,I,F]',
|
|
93 '[CR1]1[CR1][CR1][CR1][CR1][CR1][CR1]1',
|
|
94 '[CR1]1[CR1][CR1]cc[CR1][CR1]1',
|
|
95 '[CR2]1[CR2][CR2][CR2][CR2][CR2][CR2][CR2]1',
|
|
96 '[CR2]1[CR2][CR2]cc[CR2][CR2][CR2]1',
|
|
97 '[CH2R2]1N[CH2R2][CH2R2][CH2R2][CH2R2][CH2R2]1',
|
|
98 '[CH2R2]1N[CH2R2][CH2R2][CH2R2][CH2R2][CH2R2][CH2R2]1',
|
|
99 'C#C',
|
|
100 '[OR2,NR2]@[CR2]@[CR2]@[OR2,NR2]@[CR2]@[CR2]@[OR2,NR2]',
|
|
101 '[$([N+R]),$([n+R]),$([N+]=C)][O-]',
|
|
102 '[C,c]=N[OH]',
|
|
103 '[C,c]=NOC=O',
|
|
104 '[C,c](=O)[CX4,CR0X3,O][C,c](=O)',
|
|
105 'c1ccc2c(c1)ccc(=O)o2',
|
|
106 '[O+,o+,S+,s+]',
|
|
107 'N=C=O',
|
|
108 '[NX3,NX4][F,Cl,Br,I]',
|
|
109 'c1ccccc1OC(=O)[#6]',
|
|
110 '[CR0]=[CR0][CR0]=[CR0]',
|
|
111 '[C+,c+,C-,c-]',
|
|
112 'N=[N+]=[N-]',
|
|
113 'C12C(NC(N1)=O)CSC2',
|
|
114 'c1c([OH])c([OH,NH2,NH])ccc1',
|
|
115 'P',
|
|
116 '[N,O,S]C#N',
|
|
117 'C=C=O',
|
|
118 '[Si][F,Cl,Br,I]',
|
|
119 '[SX2]O',
|
|
120 '[SiR0,CR0](c1ccccc1)(c2ccccc2)(c3ccccc3)',
|
|
121 'O1CCCCC1OC2CCC3CCCCC3C2',
|
|
122 'N=[CR0][N,n,O,S]',
|
|
123 '[cR2]1[cR2][cR2]([Nv3X3,Nv4X4])[cR2][cR2][cR2]1[cR2]2[cR2][cR2][cR2]([Nv3X3,Nv4X4])[cR2][cR2]2',
|
|
124 'C=[C!r]C#N',
|
|
125 '[cR2]1[cR2]c([N+0X3R0,nX3R0])c([N+0X3R0,nX3R0])[cR2][cR2]1',
|
|
126 '[cR2]1[cR2]c([N+0X3R0,nX3R0])[cR2]c([N+0X3R0,nX3R0])[cR2]1',
|
|
127 '[cR2]1[cR2]c([N+0X3R0,nX3R0])[cR2][cR2]c1([N+0X3R0,nX3R0])',
|
|
128 '[OH]c1ccc([OH,NH2,NH])cc1',
|
|
129 'c1ccccc1OC(=O)O',
|
|
130 '[SX2H0][N]',
|
|
131 'c12ccccc1(SC(S)=N2)',
|
|
132 'c12ccccc1(SC(=S)N2)',
|
|
133 'c1nnnn1C=O',
|
|
134 's1c(S)nnc1NC=O',
|
|
135 'S1C=CSC1=S',
|
|
136 'C(=O)Onnn',
|
|
137 'OS(=O)(=O)C(F)(F)F',
|
|
138 'N#CC[OH]',
|
|
139 'N#CC(=O)',
|
|
140 'S(=O)(=O)C#N',
|
|
141 'N[CH2]C#N',
|
|
142 'C1(=O)NCC1',
|
|
143 'S(=O)(=O)[O-,OH]',
|
|
144 'NC[F,Cl,Br,I]',
|
|
145 'C=[C!r]O',
|
|
146 '[NX2+0]=[O+0]',
|
|
147 '[OR0,NR0][OR0,NR0]',
|
|
148 'C(=O)O[C,H1].C(=O)O[C,H1].C(=O)O[C,H1]',
|
|
149 '[CX2R0][NX3R0]',
|
|
150 'c1ccccc1[C;!R]=[C;!R]c2ccccc2',
|
|
151 '[NX3R0,NX4R0,OR0,SX2R0][CX4][NX3R0,NX4R0,OR0,SX2R0]',
|
|
152 '[s,S,c,C,n,N,o,O]~[n+,N+](~[s,S,c,C,n,N,o,O])(~[s,S,c,C,n,N,o,O])~[s,S,c,C,n,N,o,O]',
|
|
153 '[s,S,c,C,n,N,o,O]~[nX3+,NX3+](~[s,S,c,C,n,N])~[s,S,c,C,n,N]',
|
|
154 '[*]=[N+]=[*]',
|
|
155 '[SX3](=O)[O-,OH]',
|
|
156 'N#N',
|
|
157 'F.F.F.F',
|
|
158 '[R0;D2][R0;D2][R0;D2][R0;D2]',
|
|
159 '[cR,CR]~C(=O)NC(=O)~[cR,CR]',
|
|
160 'C=!@CC=[O,S]',
|
|
161 '[#6,#8,#16][C,c](=O)O[C,c]',
|
|
162 'c[C;R0](=[O,S])[C,c]',
|
|
163 'c[SX2][C;!R]',
|
|
164 'C=C=C',
|
|
165 'c1nc([F,Cl,Br,I,S])ncc1',
|
|
166 'c1ncnc([F,Cl,Br,I,S])c1',
|
|
167 'c1nc(c2c(n1)nc(n2)[F,Cl,Br,I])',
|
|
168 '[C,c]S(=O)(=O)c1ccc(cc1)F',
|
|
169 '[15N]',
|
|
170 '[13C]',
|
|
171 '[18O]',
|
|
172 '[34S]'
|
|
173 ]
|
|
174
|
|
175 StructuralAlerts = []
|
|
176 for smarts in StructuralAlertSmarts:
|
|
177 StructuralAlerts.append(Chem.MolFromSmarts(smarts))
|
|
178
|
|
179
|
|
180 # ADS parameters for the 8 molecular properties: [row][column]
|
|
181 # rows[8]: MW, ALOGP, HBA, HBD, PSA, ROTB, AROM, ALERTS
|
|
182 # columns[7]: A, B, C, D, E, F, DMAX
|
|
183 # ALOGP parameters from Gregory Gerebtzoff (2012, Roche)
|
|
184 pads1 = [ [2.817065973, 392.5754953, 290.7489764, 2.419764353, 49.22325677, 65.37051707, 104.9805561],
|
|
185 [0.486849448, 186.2293718, 2.066177165, 3.902720615, 1.027025453, 0.913012565, 145.4314800],
|
|
186 [2.948620388, 160.4605972, 3.615294657, 4.435986202, 0.290141953, 1.300669958, 148.7763046],
|
|
187 [1.618662227, 1010.051101, 0.985094388, 0.000000001, 0.713820843, 0.920922555, 258.1632616],
|
|
188 [1.876861559, 125.2232657, 62.90773554, 87.83366614, 12.01999824, 28.51324732, 104.5686167],
|
|
189 [0.010000000, 272.4121427, 2.558379970, 1.565547684, 1.271567166, 2.758063707, 105.4420403],
|
|
190 [3.217788970, 957.7374108, 2.274627939, 0.000000001, 1.317690384, 0.375760881, 312.3372610],
|
|
191 [0.010000000, 1199.094025, -0.09002883, 0.000000001, 0.185904477, 0.875193782, 417.7253140] ]
|
|
192 # ALOGP parameters from the original publication
|
|
193 pads2 = [ [2.817065973, 392.5754953, 290.7489764, 2.419764353, 49.22325677, 65.37051707, 104.9805561],
|
|
194 [3.172690585, 137.8624751, 2.534937431, 4.581497897, 0.822739154, 0.576295591, 131.3186604],
|
|
195 [2.948620388, 160.4605972, 3.615294657, 4.435986202, 0.290141953, 1.300669958, 148.7763046],
|
|
196 [1.618662227, 1010.051101, 0.985094388, 0.000000001, 0.713820843, 0.920922555, 258.1632616],
|
|
197 [1.876861559, 125.2232657, 62.90773554, 87.83366614, 12.01999824, 28.51324732, 104.5686167],
|
|
198 [0.010000000, 272.4121427, 2.558379970, 1.565547684, 1.271567166, 2.758063707, 105.4420403],
|
|
199 [3.217788970, 957.7374108, 2.274627939, 0.000000001, 1.317690384, 0.375760881, 312.3372610],
|
|
200 [0.010000000, 1199.094025, -0.09002883, 0.000000001, 0.185904477, 0.875193782, 417.7253140] ]
|
|
201
|
|
202 def ads(x, a, b, c, d, e, f, dmax):
|
|
203 return ((a+(b/(1+exp(-1*(x-c+d/2)/e))*(1-1/(1+exp(-1*(x-c-d/2)/f))))) / dmax)
|
|
204
|
|
205 def properties(mol):
|
|
206 """
|
|
207 Calculates the properties that are required to calculate the QED descriptor.
|
|
208 """
|
|
209 matches = []
|
|
210 if (mol is None):
|
|
211 raise WrongArgument("properties(mol)", "mol argument is \'None\'")
|
|
212 x = [0] * 9
|
|
213 x[0] = Descriptors.MolWt(mol) # MW
|
|
214 x[1] = Descriptors.MolLogP(mol) # ALOGP
|
|
215 for hba in Acceptors: # HBA
|
|
216 if (mol.HasSubstructMatch(hba)):
|
|
217 matches = mol.GetSubstructMatches(hba)
|
|
218 x[2] += len(matches)
|
|
219 x[3] = Descriptors.NumHDonors(mol) # HBD
|
|
220 x[4] = Descriptors.TPSA(mol) # PSA
|
|
221 x[5] = Descriptors.NumRotatableBonds(mol) # ROTB
|
|
222 x[6] = Chem.GetSSSR(Chem.DeleteSubstructs(deepcopy(mol), AliphaticRings)) # AROM
|
|
223 for alert in StructuralAlerts: # ALERTS
|
|
224 if (mol.HasSubstructMatch(alert)): x[7] += 1
|
|
225 ro5_failed = 0
|
|
226 if x[3] > 5:
|
|
227 ro5_failed += 1 #HBD
|
|
228 if x[2] > 10:
|
|
229 ro5_failed += 1 #HBA
|
|
230 if x[0] >= 500:
|
|
231 ro5_failed += 1
|
|
232 if x[1] > 5:
|
|
233 ro5_failed += 1
|
|
234 x[8] = ro5_failed
|
|
235 return x
|
|
236
|
|
237
|
|
238 def qed(w, p, gerebtzoff):
|
|
239 d = [0.00] * 8
|
|
240 if gerebtzoff:
|
|
241 for i in range(0, 8):
|
|
242 d[i] = ads(p[i], pads1[i][0], pads1[i][1], pads1[i][2], pads1[i][3], pads1[i][4], pads1[i][5], pads1[i][6])
|
|
243 else:
|
|
244 for i in range(0, 8):
|
|
245 d[i] = ads(p[i], pads2[i][0], pads2[i][1], pads2[i][2], pads2[i][3], pads2[i][4], pads2[i][5], pads2[i][6])
|
|
246 t = 0.0
|
|
247 for i in range(0, 8):
|
|
248 t += w[i] * log(d[i])
|
|
249 return (exp(t / sum(w)))
|
|
250
|
|
251
|
|
252 def weights_max(mol, gerebtzoff = True, props = False):
|
|
253 """
|
|
254 Calculates the QED descriptor using maximal descriptor weights.
|
|
255 If props is specified we skip the calculation step and use the props-list of properties.
|
|
256 """
|
|
257 if not props:
|
|
258 props = properties(mol)
|
|
259 return qed([0.50, 0.25, 0.00, 0.50, 0.00, 0.50, 0.25, 1.00], props, gerebtzoff)
|
|
260
|
|
261
|
|
262 def weights_mean(mol, gerebtzoff = True, props = False):
|
|
263 """
|
|
264 Calculates the QED descriptor using average descriptor weights.
|
|
265 If props is specified we skip the calculation step and use the props-list of properties.
|
|
266 """
|
|
267 if not props:
|
|
268 props = properties(mol)
|
|
269 return qed([0.66, 0.46, 0.05, 0.61, 0.06, 0.65, 0.48, 0.95], props, gerebtzoff)
|
|
270
|
|
271
|
|
272 def weights_none(mol, gerebtzoff = True, props = False):
|
|
273 """
|
|
274 Calculates the QED descriptor using unit weights.
|
|
275 If props is specified we skip the calculation step and use the props-list of properties.
|
|
276 """
|
|
277 if not props:
|
|
278 props = properties(mol)
|
|
279 return qed([1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00], props, gerebtzoff)
|
|
280
|
|
281
|
|
282 def default(mol, gerebtzoff = True):
|
|
283 """
|
|
284 Calculates the QED descriptor using average descriptor weights and Gregory Gerebtzoff parameters.
|
|
285 """
|
|
286 return weights_mean(mol, gerebtzoff)
|
|
287
|
|
288
|
|
289 if __name__ == "__main__":
|
|
290 parser = argparse.ArgumentParser()
|
|
291 parser.add_argument('-i', '--input', required=True, help='path to the input file name')
|
|
292 parser.add_argument("-m", "--method", dest="method", choices=['max', 'mean', 'unweighted'],
|
|
293 default="mean",
|
|
294 help="Specify the method you want to use.")
|
|
295
|
|
296 parser.add_argument('-o', '--outfile', type=argparse.FileType('w+'),
|
|
297 default=sys.stdout, help="path to the result file, default it sdtout")
|
|
298
|
|
299 parser.add_argument("--header", dest="header", action="store_true",
|
|
300 default=False,
|
|
301 help="Write header line.")
|
|
302
|
|
303
|
|
304 args = parser.parse_args()
|
|
305
|
|
306 # Elucidate filetype and open supplier
|
|
307 ifile = os.path.abspath(args.input)
|
|
308 if not os.path.isfile(ifile):
|
|
309 print "Error: ", ifile, " is not a file or cannot be found."
|
|
310 sys.exit(1)
|
|
311 if not os.path.exists(ifile):
|
|
312 print "Error: ", ifile, " does not exist or cannot be found."
|
|
313 sys.exit(1)
|
|
314 if not os.access(ifile, os.R_OK):
|
|
315 print "Error: ", ifile, " is not readable."
|
|
316 sys.exit(1)
|
|
317
|
|
318 filetype = check_filetype(ifile)
|
|
319
|
|
320
|
|
321 """
|
|
322 We want to store the original SMILES in the output. So in case of a SMILES file iterate over the file and convert each line separate.
|
|
323 """
|
|
324
|
|
325 if filetype == 'sdf':
|
|
326 supplier = Chem.SDMolSupplier(ifile)
|
|
327 # Process file
|
|
328 if args.header:
|
|
329 args.outfile.write("MW\tALOGP\tHBA\tHBD\tPSA\tROTB\tAROM\tALERTS\tLRo5\tQED\tNAME\n")
|
|
330 count = 0
|
|
331 for mol in supplier:
|
|
332 count += 1
|
|
333 if mol is None:
|
|
334 print "Warning: skipping molecule ", count, " and continuing with next."
|
|
335 continue
|
|
336 props = properties(mol)
|
|
337
|
|
338 if args.method == 'max':
|
|
339 calc_qed = weights_max(mol, True, props)
|
|
340 elif args.method == 'unweighted':
|
|
341 calc_qed = weights_none(mol, True, props)
|
|
342 else:
|
|
343 calc_qed = weights_mean(mol, True, props)
|
|
344
|
|
345 args.outfile.write( "%.2f\t%.3f\t%d\t%d\t%.2f\t%d\t%d\t%d\t%s\t%.3f\t%-s\n" % (
|
|
346 props[0],
|
|
347 props[1],
|
|
348 props[2],
|
|
349 props[3],
|
|
350 props[4],
|
|
351 props[5],
|
|
352 props[6],
|
|
353 props[7],
|
|
354 props[8],
|
|
355 calc_qed,
|
|
356 mol.GetProp("_Name"),
|
|
357 ))
|
|
358 elif filetype == 'smi':
|
|
359 supplier = Chem.SmilesMolSupplier(ifile, " \t", 0, 1, False, True)
|
|
360
|
|
361 # Process file
|
|
362 if args.header:
|
|
363 args.outfile.write("MW\tALOGP\tHBA\tHBD\tPSA\tROTB\tAROM\tALERTS\tLRo5\tQED\tNAME\tSMILES\n")
|
|
364 count = 0
|
|
365 for line in open(ifile):
|
|
366 tokens = line.strip().split('\t')
|
|
367 if len(tokens) > 1:
|
|
368 smiles, title = tokens
|
|
369 else:
|
|
370 smiles = tokens[0]
|
|
371 title = ''
|
|
372 mol = Chem.MolFromSmiles(smiles)
|
|
373 count += 1
|
|
374 if mol is None:
|
|
375 print "Warning: skipping molecule ", count, " and continuing with next."
|
|
376 continue
|
|
377 props = properties(mol)
|
|
378
|
|
379 if args.method == 'max':
|
|
380 calc_qed = weights_max(mol, True, props)
|
|
381 elif args.method == 'unweighted':
|
|
382 calc_qed = weights_none(mol, True, props)
|
|
383 else:
|
|
384 calc_qed = weights_mean(mol, True, props)
|
|
385
|
|
386 args.outfile.write( "%.2f\t%.3f\t%d\t%d\t%.2f\t%d\t%d\t%d\t%s\t%.3f\t%-s\t%s\n" % (
|
|
387 props[0],
|
|
388 props[1],
|
|
389 props[2],
|
|
390 props[3],
|
|
391 props[4],
|
|
392 props[5],
|
|
393 props[6],
|
|
394 props[7],
|
|
395 props[8],
|
|
396 calc_qed,
|
|
397 title,
|
|
398 smiles
|
|
399 ))
|
|
400
|
|
401 else:
|
|
402 print "Error: unknown file extension: ", extension
|
|
403 sys.exit(1)
|
|
404
|
|
405 sys.exit(0)
|