1
2
3
4
5
6 from __future__ import print_function
7 from rdkit import RDLogger
8 logger = RDLogger.logger()
9 from rdkit import Chem,Geometry
10 import numpy
11 from rdkit.Numerics import Alignment
12 from rdkit.Chem.Subshape import SubshapeObjects
13
22
24 for i in range(len(pts)):
25 if orderedTraversal:
26 jStart=i+1
27 else:
28 jStart=0
29 for j in range(jStart,len(pts)):
30 if j==i:
31 continue
32 if orderedTraversal:
33 kStart=j+1
34 else:
35 kStart=0
36 for k in range(j+1,len(pts)):
37 if k==i or k==j:
38 continue
39 yield (i,j,k)
40
44
45
57
58
63 from rdkit.ML.Cluster import Butina
64 dists = []
65 for i in range(len(alignments)):
66 TransformMol(mol,alignments[i].transform,newConfId=tempConfId)
67 shapeI=builder.GenerateSubshapeShape(mol,tempConfId,addSkeleton=False)
68 for j in range(i):
69 TransformMol(mol,alignments[j].transform,newConfId=tempConfId+1)
70 shapeJ=builder.GenerateSubshapeShape(mol,tempConfId+1,addSkeleton=False)
71 d = GetShapeShapeDistance(shapeI,shapeJ,distMetric)
72 dists.append(d)
73 mol.RemoveConformer(tempConfId+1)
74 mol.RemoveConformer(tempConfId)
75 clusts=Butina.ClusterData(dists,len(alignments),neighborTol,isDistData=True)
76 res = [alignments[x[0]] for x in clusts]
77 return res
78
95
97 triangleRMSTol=1.0
98 distMetric=SubshapeDistanceMetric.PROTRUDE
99 shapeDistTol=0.2
100 numFeatThresh=3
101 dirThresh=2.6
102 edgeTol=6.0
103
104
105 coarseGridToleranceMult=1.0
106 medGridToleranceMult=1.0
107
109 """ this is a generator function returning the possible triangle
110 matches between the two shapes
111 """
112 ssdTol = (self.triangleRMSTol**2)*9
113 res = []
114 tgtPts = target.skelPts
115 queryPts = query.skelPts
116 tgtLs = {}
117 for i in range(len(tgtPts)):
118 for j in range(i+1,len(tgtPts)):
119 l2 = (tgtPts[i].location-tgtPts[j].location).LengthSq()
120 tgtLs[(i,j)]=l2
121 queryLs = {}
122 for i in range(len(queryPts)):
123 for j in range(i+1,len(queryPts)):
124 l2 = (queryPts[i].location-queryPts[j].location).LengthSq()
125 queryLs[(i,j)]=l2
126 compatEdges={}
127 tol2 = self.edgeTol*self.edgeTol
128 for tk,tv in tgtLs.iteritems():
129 for qk,qv in queryLs.iteritems():
130 if abs(tv-qv)<tol2:
131 compatEdges[(tk,qk)]=1
132 seqNo=0
133 for tgtTri in _getAllTriangles(tgtPts,orderedTraversal=True):
134 tgtLocs=[tgtPts[x].location for x in tgtTri]
135 for queryTri in _getAllTriangles(queryPts,orderedTraversal=False):
136 if compatEdges.has_key(((tgtTri[0],tgtTri[1]),(queryTri[0],queryTri[1]))) and \
137 compatEdges.has_key(((tgtTri[0],tgtTri[2]),(queryTri[0],queryTri[2]))) and \
138 compatEdges.has_key(((tgtTri[1],tgtTri[2]),(queryTri[1],queryTri[2]))):
139 queryLocs=[queryPts[x].location for x in queryTri]
140 ssd,tf = Alignment.GetAlignmentTransform(tgtLocs,queryLocs)
141 if ssd<=ssdTol:
142 alg = SubshapeAlignment()
143 alg.transform=tf
144 alg.triangleSSD=ssd
145 alg.targetTri=tgtTri
146 alg.queryTri=queryTri
147 alg._seqNo=seqNo
148 seqNo+=1
149 yield alg
150
152 nMatched=0
153 for i in range(3):
154 tgtFeats = targetPts[alignment.targetTri[i]].molFeatures
155 qFeats = queryPts[alignment.queryTri[i]].molFeatures
156 if not tgtFeats and not qFeats:
157 nMatched+=1
158 else:
159 for j,jFeat in enumerate(tgtFeats):
160 if jFeat in qFeats:
161 nMatched+=1
162 break
163 if nMatched>=self.numFeatThresh:
164 break
165 return nMatched>=self.numFeatThresh
166
168 i = 0
169 targetPts = target.skelPts
170 queryPts = query.skelPts
171 while i<len(alignments):
172 alg = alignments[i]
173 if not self._checkMatchFeatures(targetPts,queryPts,alg):
174 if pruneStats is not None:
175 pruneStats['features']=pruneStats.get('features',0)+1
176 del alignments[i]
177 else:
178 i+=1
179
197
199 i = 0
200 tgtPts = target.skelPts
201 queryPts = query.skelPts
202 while i<len(alignments):
203 if not self._checkMatchDirections(tgtPts,queryPts,alignments[i]):
204 if pruneStats is not None:
205 pruneStats['direction']=pruneStats.get('direction',0)+1
206 del alignments[i]
207 else:
208 i+=1
209
221
222 - def _checkMatchShape(self,targetMol,target,queryMol,query,alignment,builder,
223 targetConf,queryConf,pruneStats=None,tConfId=1001):
224 matchOk=True
225 TransformMol(queryMol,alignment.transform,confId=queryConf,newConfId=tConfId)
226 oSpace=builder.gridSpacing
227 builder.gridSpacing=oSpace*2
228 coarseGrid=builder.GenerateSubshapeShape(queryMol,tConfId,addSkeleton=False)
229 d = GetShapeShapeDistance(coarseGrid,target.coarseGrid,self.distMetric)
230 if d>self.shapeDistTol*self.coarseGridToleranceMult:
231 matchOk=False
232 if pruneStats is not None:
233 pruneStats['coarseGrid']=pruneStats.get('coarseGrid',0)+1
234 else:
235 builder.gridSpacing=oSpace*1.5
236 medGrid=builder.GenerateSubshapeShape(queryMol,tConfId,addSkeleton=False)
237 d = GetShapeShapeDistance(medGrid,target.medGrid,self.distMetric)
238 if d>self.shapeDistTol*self.medGridToleranceMult:
239 matchOk=False
240 if pruneStats is not None:
241 pruneStats['medGrid']=pruneStats.get('medGrid',0)+1
242 else:
243 builder.gridSpacing=oSpace
244 fineGrid=builder.GenerateSubshapeShape(queryMol,tConfId,addSkeleton=False)
245 d = GetShapeShapeDistance(fineGrid,target,self.distMetric)
246
247 if d>self.shapeDistTol:
248 matchOk=False
249 if pruneStats is not None:
250 pruneStats['fineGrid']=pruneStats.get('fineGrid',0)+1
251 alignment.shapeDist=d
252 queryMol.RemoveConformer(tConfId)
253 builder.gridSpacing=oSpace
254 return matchOk
255
256 - def PruneMatchesUsingShape(self,targetMol,target,queryMol,query,builder,
257 alignments,tgtConf=-1,queryConf=-1,
258 pruneStats=None):
259 if not hasattr(target,'medGrid'):
260 self._addCoarseAndMediumGrids(targetMol,target,tgtConf,builder)
261
262 logger.info("Shape-based Pruning")
263 i=0
264 nOrig = len(alignments)
265 nDone=0
266 while i < len(alignments):
267 removeIt=False
268 alg = alignments[i]
269 nDone+=1
270 if not nDone%100:
271 nLeft = len(alignments)
272 logger.info(' processed %d of %d. %d alignments remain'%((nDone,
273 nOrig,
274 nLeft)))
275 if not self._checkMatchShape(targetMol,target,queryMol,query,alg,builder,
276 targetConf=tgtConf,queryConf=queryConf,
277 pruneStats=pruneStats):
278 del alignments[i]
279 else:
280 i+=1
281
282 - def GetSubshapeAlignments(self,targetMol,target,queryMol,query,builder,
283 tgtConf=-1,queryConf=-1,pruneStats=None):
284 import time
285 if pruneStats is None:
286 pruneStats={}
287 logger.info("Generating triangle matches")
288 t1=time.time()
289 res = [x for x in self.GetTriangleMatches(target,query)]
290 t2=time.time()
291 logger.info("Got %d possible alignments in %.1f seconds"%(len(res),t2-t1))
292 pruneStats['gtm_time']=t2-t1
293 if builder.featFactory:
294 logger.info("Doing feature pruning")
295 t1 = time.time()
296 self.PruneMatchesUsingFeatures(target,query,res,pruneStats=pruneStats)
297 t2 = time.time()
298 pruneStats['feats_time']=t2-t1
299 logger.info("%d possible alignments remain. (%.1f seconds required)"%(len(res),t2-t1))
300 logger.info("Doing direction pruning")
301 t1 = time.time()
302 self.PruneMatchesUsingDirection(target,query,res,pruneStats=pruneStats)
303 t2 = time.time()
304 pruneStats['direction_time']=t2-t1
305 logger.info("%d possible alignments remain. (%.1f seconds required)"%(len(res),t2-t1))
306 t1 = time.time()
307 self.PruneMatchesUsingShape(targetMol,target,queryMol,query,builder,res,
308 tgtConf=tgtConf,queryConf=queryConf,
309 pruneStats=pruneStats)
310 t2 = time.time()
311 pruneStats['shape_time']=t2-t1
312 return res
313
314 - def __call__(self,targetMol,target,queryMol,query,builder,
315 tgtConf=-1,queryConf=-1,pruneStats=None):
316 for alignment in self.GetTriangleMatches(target,query):
317 if builder.featFactory and \
318 not self._checkMatchFeatures(target.skelPts,query.skelPts,alignment):
319 if pruneStats is not None:
320 pruneStats['features']=pruneStats.get('features',0)+1
321 continue
322 if not self._checkMatchDirections(target.skelPts,query.skelPts,alignment):
323 if pruneStats is not None:
324 pruneStats['direction']=pruneStats.get('direction',0)+1
325 continue
326
327 if not hasattr(target,'medGrid'):
328 self._addCoarseAndMediumGrids(targetMol,target,tgtConf,builder)
329
330 if not self._checkMatchShape(targetMol,target,queryMol,query,alignment,builder,
331 targetConf=tgtConf,queryConf=queryConf,
332 pruneStats=pruneStats):
333 continue
334
335 yield alignment
336
337
338 if __name__=='__main__':
339 from rdkit.six.moves import cPickle
340 tgtMol,tgtShape = cPickle.load(file('target.pkl','rb'))
341 queryMol,queryShape = cPickle.load(file('query.pkl','rb'))
342 builder = cPickle.load(file('builder.pkl','rb'))
343 aligner = SubshapeAligner()
344 algs = aligner.GetSubshapeAlignments(tgtMol,tgtShape,queryMol,queryShape,builder)
345 print(len(algs))
346
347 from rdkit.Chem.PyMol import MolViewer
348 v = MolViewer()
349 v.ShowMol(tgtMol,name='Target',showOnly=True)
350 v.ShowMol(queryMol,name='Query',showOnly=False)
351 SubshapeObjects.DisplaySubshape(v,tgtShape,'target_shape',color=(.8,.2,.2))
352 SubshapeObjects.DisplaySubshape(v,queryShape,'query_shape',color=(.2,.2,.8))
353