source: orange-bioinformatics/orangecontrib/bio/obiChem.py @ 1873:0810c5708cc5

Revision 1873:0810c5708cc5, 49.7 KB checked in by Ales Erjavec <ales.erjavec@…>, 6 months ago (diff)

Moved '_bioinformatics' into orangecontrib namespace.

RevLine 
[597]1"""A library for searching frequent molecular fragments (substructures) based on the
2Mining Molecular Fragments: Finding Relevant Substructures of Molecules
3Christian Borgelt and Michael R. Berthold.
4Classes (see their corresponding __doc__ strings for further detail):
5    Fragment        : Representation of the fragment
6    FragmentMiner   : The main class that does the search
7    Fragmenter      : A class that is used to fragment an ExampleTable
8    FragmentBasedLearner    : A learner wrapper class that first runs the molecular fragmentation on the data
9"""
[949]10try:
11    from openbabel import OBMol, OBAtom, OBBond, OBSmartsPattern, OBConversion, OBMolAtomIter, OBMolBondIter, OBAtomBondIter
12except ImportError:
13    pass
[597]14##from pybel import *
15from copy import deepcopy
16import orange
[606]17import sys, os
[597]18
19
20debug=False
21
22class Atom(object):
23    def __init__(self, molecule, atomicNum=6, aromaticFlag=False, atomIndex=0):
24        self.molecule=molecule
25        self.atomicNum=atomicNum
26        self.aromaticFlag=aromaticFlag
27        self.atomIndex=atomIndex
28        self.extendedIndex=0
29        self.lastExtendedBondOrder=1
30        self.lastExtendedAtomicNum=1
31        self.molecule.AddAtom(self)
32    def GetAtomicNum(self): return self.atomicNum
33    def IsAromatic(self): return self.aromaticFlag
34    def Match(self, atom):
35        #Match this atom to an OBAtom
36        if self.atomicNum==atom.GetAtomicNum() and self.IsAromatic()==atom.IsAromatic():
37            return True
38        else:
39            return False
40    def GetIdx(self):
41        return self.atomIndex
42
43class Bond(object):
44    def __init__(self, molecule, atom1=None, atom2=None, bondOrder=1):
45        self.molecule=molecule
46        self.atom1=atom1
47        self.atom2=atom2
48        self.bondOrder=bondOrder
49        self.molecule.AddBond(self)
50    def GetBondOrder(self): return self.bondOrder
51    def GetNbrAtom(self, atom):
52        if atom==self.atom1:
53            return self.atom2
54        elif atom==self. atom2:
55            return  self.atom1
56        else:
57            raise Exception("Atom does not belong to this bond")
58    def Match(self, bond):
59        #Match this bond to an OBBond
60        if self.GetbondOrder()==bond.GetBondOrder():
61            return True
62        else:
63            return False
64
65class Ring(object):
66    def __init__(self, molecule, obRing, embeding):
67        self.molecule=molecule
68        molecule.AddRing(self)
69        self.ringAtoms=[]
70        self.isAromatic=obRing.IsAromatic()
71        reverseEmbedingDict=embeding.GetReverseDict()
72        newAtoms={}
73        self.extendedIndex=oldNextExtendedIndex=molecule.nextExtendedIndex
74        for ind1, ind2 in GetRingPairs(obRing._path):
75            atom1=embeding.GetReverseEmbededAtom(ind1, None) or newAtoms.get(ind1, None)
76            if not atom1:
77                atom1=Atom(molecule, embeding.molecule.GetAtom(ind1).GetAtomicNum(), obRing.IsAromatic())
78                newAtoms[ind1]=atom1
79                embeding[atom1.atomIndex]=ind1
80            atom2=embeding.GetReverseEmbededAtom(ind2, None) or newAtoms.get(ind2, None)
81            if not atom2:
82                atom2=Atom(molecule, embeding.molecule.GetAtom(ind2).GetAtomicNum(), obRing.IsAromatic())
83                newAtoms[ind2]=atom2
84                embeding[atom2.atomIndex]=ind2
85            Bond(molecule, atom1, atom2, embeding.molecule.GetBond(ind1, ind2).GetBondOrder())
86            self.ringAtoms.append(atom1)
87        for atom in newAtoms.values():
88            atom.extendedIndex=oldNextExtendedIndex
89        self.extendedIndex=oldNextExtendedIndex
90        molecule.nextExtendedIndex=oldNextExtendedIndex+1
91    def Size(self):
92        return len(self.ringAtoms)
93    def __getattr__(self, name):
94        if name == "_path":
95            return [a.GetIdx() for a in self.ringAtoms]
96        raise AttributeError 
97
98class Molecule(object):
99    def __init__(self):
100        self.atoms=[]
101        self.bonds=[]
102        self.rings=[]
103        self.nextExtendedIndex=0
104
105    def __deepcopy__(self, memo):
106        mol=Molecule()
107        memo[id(self)]=mol
108        mol.__in_place_deepcopy__(self, memo)
109        return mol
110    def __in_place_deepcopy__(self, mol, memo):
111        if debug: print "Copying molecule"
112        self.atoms=deepcopy(mol.atoms, memo)
113        #memo[mol.atoms]=self.atoms
114        self.bonds=deepcopy(mol.bonds, memo)
115        #memo[mol.bonds]=self.bonds
116        self.rings=deepcopy(mol.rings, memo)
117        self.nextExtendedIndex=mol.nextExtendedIndex
118        #memo[mol.rings]=self.rings"""
119       
120    def AddAtom(self, atom):
121        if atom in self.atoms:
122            raise Exception("Atom already present")
123        self.atoms.append(atom)
124        atom.atomIndex=len(self.atoms)-1
125        atom.extendedIndex=self.nextExtendedIndex
126        self.nextExtendedIndex+=1
127    def AddBond(self, bond):
128        if bond in self.bonds:
129            raise Exception("Bond already present")
130        self.bonds.append(bond)
131        bond.bondIndex=len(self.bonds)-1
132    def AddRing(self, ring):
133        if ring in self.rings:
134            raise Exception("Ring already present")
135        self.rings.append(ring)
136        ring.ringIndex=len(self.rings)-1
137    def GetBond(self, atom1, atom2):
138        for bond in self.bonds:
139            if bond.atom1==atom1 and bond.atom2==atom2:
140                return bond
141            if bond.atom1==atom2 and bond.atom2==atom1:
142                return bond
143        return None
144    def GetAtom(self, index):
145        return self.atoms[index]
146   
147def cmpAtoms(atom1, atom2):
148    return atom1.GetAtomicNum()==atom2.GetAtomicNum() and atom1.IsAromatic()==atom2.IsAromatic()
149def cmpAtomBonds(bond1, bond2):
150    return bond1.IsAromatic() and bond2.IsAromatic() or bond1.GetBondOrder()==bond2.GetBondOrder()
151
152class Embeding(dict):
153    def __init__(self, embeding={}, molecule=None, fragment=None):
154        dict.__init__(self, embeding)
155        self.molecule=molecule
156        self.fragment=fragment
157        if embeding.__class__==Embeding:
158            if not self.molecule:
159                self.molecule=embeding.molecule
160            if not self.fragment:
161                self.fragment=embeding.fragment
162    def __deepcopy__(self, memo):
163        if id(self.fragment) in memo:
164            return Embeding(self, fragment=memo[id(self.fragment)])
165        else:
166            return Embeding(self)
167    def GetEmbededAtom(self, atomInd):
168        return self.molecule.GetAtom(self.__getitem__(atomInd))
169    def GetReverseEmbededAtom(self, atomInd, default=None):
170        rev=self.GetReverseDict()
171        if atomInd in rev:
172            return self.fragment.atoms[rev[atomInd]]
173        else:
174            return default
175    def GetReverseDict(self):
176        return dict([(value, key) for key, value in self.items()])
177
178def GetRingPairs(list):
179    return reduce(lambda a,b:a+[(a[-1][1],b)], list[1:], [(list[-1],list[0])])
180
181class FragmentExtension(object):
182    def __init__(self, startAtomInd, embedings):
183        self.startAtomInd=startAtomInd
184        self.embedings=embedings
185    def IsEquivalent(self, extension):
186        pass
187    def MergeFrom(self, extension):
188        pass
189    def Extend(self, fragment):
190        pass
191
192class FragmentExtensionByAtom(FragmentExtension):
193    def __init__(self, atomicNum, aromatic, atomIndices, embedings):
194        self.atomicNum=atomicNum
195        self.aromatic=aromatic
196        self.atomIndices=atomIndices
197        self.embedings=embedings
198    def IsEquivalent(self, extension):
199        return self.atomicNum==extension.atomicNum and self.aromatic==extension.aromatic
200    def MergeFrom(self, extension):
201        self.embedings.extend(extension.embedings)
202        self.atomIndices.extend(extension.atomIndices)
203    def Extend(self, fragment):
204        atom=Atom(fragment, self.atomicNum, self.aromatic)
205        embedings=[Embeding(e, fragment=fragment)for e in self.embedings]
206        for embeding, atomInd in zip(embedings, self.atomIndices):
207            embeding[atom.atomIndex]=atomInd
208        fragment.embedings=embedings
209       
210class FragmentExtensionByBondAtom(FragmentExtension):
211    def __init__(self, startAtomInd, bondOrder, endAtoms, embedings):
212        FragmentExtension.__init__(self, startAtomInd, embedings)
213        self.bondOrder=bondOrder
214        self.endAtoms=endAtoms
215    def IsEquivalent(self, extension):
216        return self.startAtomInd==extension.startAtomInd and self.bondOrder==extension.bondOrder and cmpAtoms(self.endAtoms[0], extension.endAtoms[0])
217    def MergeFrom(self, extension):
218        self.endAtoms.extend(extension.endAtoms)
219        self.embedings.extend(extension.embedings)
220    def Extend(self, fragment):
221        if debug: print "Extending by bond and atom"
222        sAtom=fragment.atoms[self.startAtomInd]
223        atom=Atom(fragment, self.endAtoms[0].GetAtomicNum(),self.endAtoms[0].IsAromatic())
224        #atom.SetAtomicNum(self.endAtoms[0].GetAtomicNum())
225       
226        Bond(fragment, sAtom, atom, self.bondOrder)
227        sAtom.lastExtendedBondOrder=self.bondOrder
228        sAtom.lastExtendedAtomicNum=atom.GetAtomicNum()
229        fragment.lastExtendedAtomIndex=sAtom.extendedIndex
230        embedings=[Embeding(e, fragment=fragment)for e in self.embedings]
231        for embeding, endAtom in zip(embedings, self.endAtoms):
232            embeding[atom.atomIndex]=endAtom.GetIdx()
233        fragment.embedings=embedings
234       
235class FragmentExtensionByBond(FragmentExtension):
236    def __init__(self, startAtomInd, bondOrder, endAtomInd, embedings):
237        FragmentExtension.__init__(self, startAtomInd, embedings)
238        self.bondOrder=bondOrder
239        self.endAtomInd=endAtomInd
240    def IsEquivalent(self, extension):
241        return self.startAtomInd==extension.startAtomInd and self.bondOrder==extension.bondOrder and self.endAtomInd==extension.endAtomInd
242    def MergeFrom(self, extension):
243        self.embedings.extend(extension.embedings)
244    def Extend(self, fragment):
245        if debug: print "Extending by bond"
246        sAtom=fragment.atoms[self.startAtomInd]
247        eAtom=fragment.atoms[self.endAtomInd]
248        Bond(fragment, sAtom, eAtom, self.bondOrder)
249        sAtom.lastExtendedBondOrder=self.bondOrder
250        sAtom.lastExtendedAtomicNum=eAtom.GetAtomicNum()
251        fragment.lastExtendedAtomIndex=sAtom.extendedIndex
252        fragment.embedings=[Embeding(e, fragment=fragment) for e in self.embedings]
253
254class FragmentExtensionByBondRing(FragmentExtension):
255    def __init__(self, startAtomInd, bondOrder, endAtoms, rings, embedings):
256        FragmentExtension.__init__(self, startAtomInd, embedings)
257        self.bondOrder=bondOrder
258        self.endAtoms=endAtoms
259        self.rings=rings
260    def GetRingMapping(self, endAtom1, ring1, embeding1, endAtom2, ring2, embeding2):
261        rDict1=embeding1.GetReverseDict()
262        rDict2=embeding2.GetReverseDict()
263        ringMappings=[]
264        i=list(ring1._path).index(endAtom1.GetIdx())
265        c0=list(ring1._path)[i:]+list(ring1._path)[:i]
266        i=list(ring2._path).index(endAtom2.GetIdx())
267        c1=list(ring2._path)[i:]+list(ring2._path)[:i]
268        l=list(c1[1:]) #l=list(c1[:-1])
269        l.reverse()
270        c2=[c1[0]]+l #c2=[c1[-1]]+l
271        for c in [c1,c2]:
272            mapping=[]
273            for (ra11, ra12), (ra21, ra22) in zip(GetRingPairs(map(lambda i:ring1.GetParent().GetAtom(i), c0)),GetRingPairs(map(lambda i:ring2.GetParent().GetAtom(i), c))):
274                if cmpAtoms(ra12, ra22) and cmpAtomBonds(ring1.GetParent().GetBond(ra11, ra12), ring2.GetParent().GetBond(ra21, ra22)):
275                    mapping.append((ra12.GetIdx(), ra22.GetIdx()))
276                else:
277                    break
278            else:
279                ringMappings.append(dict(mapping))
280        return ringMappings
281    def IsEquivalent(self, extension):
282        if self.startAtomInd==extension.startAtomInd and self.bondOrder==extension.bondOrder and len(self.rings[0]._path)==len(extension.rings[0]._path) \
283           and self.rings[0].fingerprint==extension.rings[0].fingerprint:
284            return bool(self.GetRingMapping(self.endAtoms[0], self.rings[0], self.embedings[0], extension.endAtoms[0], extension.rings[0], extension.embedings[0]))
285        return False
286    def MergeFrom(self, extension):
287        self.endAtoms.extend(extension.endAtoms)
288        self.embedings.extend(extension.embedings)
289        self.rings.extend(extension.rings)
290    def Extend(self, fragment):
291        if debug: print "Extending by ring and bond"
292        sAtom=fragment.atoms[self.startAtomInd]
293        endAtom1=self.endAtoms[0]
294        ring1=self.rings[0]    #One ring to rule them all
295        embeding1=Embeding(self.embedings[0], fragment=fragment)
296        #Add all the atoms and bonds in the ring
297        newRing=Ring(fragment, ring1, embeding1) #Changes embeding
298        Bond(fragment, fragment.atoms[self.startAtomInd], fragment.atoms[embeding1.GetReverseDict()[endAtom1.GetIdx()]], self.bondOrder)
299        sAtom.lastExtendedAtomicNum=endAtom1.GetAtomicNum()
300        sAtom.lastExtendedBondOrder=self.bondOrder
301        fragment.lastExtendedAtomIndex=sAtom.extendedIndex
302        #Update embedings
303        #embeding1=self.embedings[0]
304        fragment.embedings=[]
305        for endAtom, ring, embeding in zip(self.endAtoms, self.rings, self.embedings):
306            #embeding[sAtom.atomIndex]=endAtom.GetIdx()
307            mappings=self.GetRingMapping(endAtom1, ring1, embeding1, endAtom, ring, embeding)
308            for map in mappings:
309                embeding=Embeding(embeding, fragment=fragment)
310                for atom in newRing.ringAtoms:
311                    embeding[atom.atomIndex]=map[embeding1[atom.atomIndex]]
312                embeding[newRing]=ring
313                fragment.embedings.append(embeding)
314       
315class FragmentExtensionByRing(FragmentExtension):
316    def __init__(self, rings, embedings):
317        FragmentExtension.__init__(self, 0, embedings)
318        self.rings=rings
319    def GetRingMapping(self, ring1, embeding1, ring2, embeding2):
320        rDict1=embeding1.GetReverseDict()
321        rDict2=embeding2.GetReverseDict()
322        ringMappings=[]
323        for i in range(ring2.Size()):
324            c1=list(ring2._path)[i:]+list(ring2._path)[:i]
325            l=list(c1[1:])
326            l.reverse()
327            c2=[c1[0]]+l
328            for c in [c1,c2]:
329                mapping=[]
330                for (ra11, ra12), (ra21, ra22) in zip(GetRingPairs(map(lambda i:ring1.GetParent().GetAtom(i), ring1._path)),GetRingPairs(map(lambda i:ring2.GetParent().GetAtom(i), c))):
331                    if cmpAtoms(ra12, ra22) and cmpAtomBonds(ring1.GetParent().GetBond(ra11, ra12), ring2.GetParent().GetBond(ra21, ra22)) and rDict1.get(ra12.GetIdx(), None)==rDict2.get(ra22.GetIdx(), None):
332                        mapping.append((ra12.GetIdx(), ra22.GetIdx()))
333                    else:
334                        break
335                else:
336                    ringMappings.append(dict(mapping))
337        return ringMappings
338    def IsEquivalent(self, extension):
339        if len(self.rings[0]._path)==len(extension.rings[0]._path) and self.rings[0].fingerprint==extension.rings[0].fingerprint:
340            return bool(self.GetRingMapping(self.rings[0], self.embedings[0], extension.rings[0], extension.embedings[0]))
341        return False
342    def MergeFrom(self, extension):
343        self.rings.extend(extension.rings)
344        self.embedings.extend(extension.embedings)
345    def Extend(self, fragment):
346        if debug: print "Extending by ring"
347        ring1=self.rings[0]
348        tmpEmbeding=Embeding(self.embedings[0],fragment=fragment)
349        newRing=Ring(fragment, ring1, tmpEmbeding) #Changes embeding
350        embeding1=self.embedings[0]
351##        fragment.lastExtendedAtomIndex=sAtom.extendedIndex
352        fragment.embedings=[]
353        for ring, embeding in zip(self.rings, self.embedings):
354            mappings=self.GetRingMapping(ring1, embeding1, ring, embeding)
355            for map in mappings:
356                embeding=Embeding(embeding, fragment=fragment)
357                for atom in newRing.ringAtoms:
358                    embeding[atom.atomIndex]=map[tmpEmbeding[atom.atomIndex]]
359                embeding[newRing]=ring
360                fragment.embedings.append(embeding)
361       
362class Fragment(Molecule):
363    """A class representing a molecular fragment
364    Methods:
365        ToOBMol()   : Returns an openbabel.OBMol object representation
366        ToSmiles()  : Returns a SMILES code representation
367        ToCanonicalSmiles() : Returns a canonical SMILES code representation
368        Support()   : Returns the support of the fragment in the active set
369        OcurrencesIn(smiles): Returns the number of times a fragment is containd
370                    in the molecule represented by the smiles code argument
371        ContainedIn(smiles) : Returns True if the fragment is present in the molecule
372                    represented by the smiles code argument
373    """
374    def __init__(self, miner=None, excludeAtomList=[]):
375        Molecule.__init__(self)
376        self.embedings=[]
377        self.miner=miner
378        self.excludeAtomList=excludeAtomList
379        self.lastExtendedAtomicNum=0
380        self.lastExtendedAtomIndex=0
[953]381        self.writer=OBConversion()
382        self.writer.SetInAndOutFormats("smi","smi")
[597]383
384    def __deepcopy__(self, memo):
385        f=Fragment()
386        memo[id(self)]=f
387        Molecule.__in_place_deepcopy__(f, self, memo)
388        f.embedings=[Embeding(e,fragment=f) for e in self.embedings]
389        f.miner=self.miner
390        f.excludeAtomList=self.excludeAtomList
391        f.lastExtendedAtomIndex=self.lastExtendedAtomIndex
392        return f
393       
394    def InitializeFragment(self, atomicNum):
395        mol=OBMol()
396        atom=Atom(self, atomicNum)
397        for mol in self.miner.GetAllMolecules():
398            for a in OBMolAtomIter(mol):
399                if atom.Match(a):
400                    self.embedings.append(Embeding({atom.atomIndex : a.GetIdx()}, molecule=mol, fragment=self))
401                   
402    def IsAtomExcluded(self, atom):
403        return atom.GetAtomicNum() in self.excludeAtomList
404   
405    def GetCandidateAtoms(self):
406        return filter(lambda a:a.extendedIndex>=self.lastExtendedAtomIndex, self.atoms)
407
408    def GetCandidateRings(self):
409        return filter(lambda r:r.extendedIndex>=self.lastExtendedAtomIndex, self.rings)
410
411    def FilterRings(self, candidateRings, embeding):
412        for ring1 in candidateRings:
413            for ring2 in filter(lambda r:type(r)==type(ring1), embeding.values()):
414                if ring1.this==ring2.this:
415                    break
416            else:
417                yield ring1
418       
419    def GetCandidatesFromAtom(self, atomInd, embeding):
420        candidates=[]
421        reverseEmbedingDict=embeding.GetReverseDict()
422        atom=embeding.GetEmbededAtom(atomInd)
423        if atom.GetParent().NumHvyAtoms()==1:
424            return candidates
425        for bond in OBAtomBondIter(atom):   #crashes if the atom is alone in a molecule
426            nbrAtom=bond.GetNbrAtom(atom)
427            if self.IsAtomExcluded(nbrAtom):
428                    continue
429            if bond.GetBondOrder()>self.atoms[atomInd].lastExtendedBondOrder or (bond.GetBondOrder()==self.atoms[atomInd].lastExtendedBondOrder and nbrAtom.GetAtomicNum()>=self.atoms[atomInd].lastExtendedAtomicNum):
430                if self.miner.addWholeRings: #Whole rings are added at the same time (no new bond can connect to an atom already in the embeding)
431                    if nbrAtom.IsInRing():
432                        if not bond.IsInRing(): #ring to ring extensions are handled in GetRingCandidatesFromRing
433                            for ring in self.FilterRings(self.miner.rings[embeding.molecule], embeding):
434                                if ring.IsMember(nbrAtom):
435                                    candidates.append(FragmentExtensionByBondRing(atomInd, bond.GetBondOrder(), [nbrAtom], [ring], [embeding]))
436##                        elif len(self.atoms)==1: #
437##                            for ring in self.FilterRings(self.miner.rings[embeding.molecule], embeding):
438##                                if ring.IsMember(nbrAtom) and ring.IsMember(atom) and ring not in addedRings:
439##                                    candidates.append(FragmentExtensionByRing([ring], [embeding]))
440##                                    addedRings.append(ring)
441                           
442                    elif nbrAtom.GetIdx() not in reverseEmbedingDict:
443                        candidates.append(FragmentExtensionByBondAtom(atomInd, bond.GetBondOrder(), [nbrAtom], [embeding]))                               
444                else: 
445                    if nbrAtom.GetIdx() in reverseEmbedingDict:
446                        if not self.GetBond(self.GetAtom(atomInd), self.GetAtom(reverseEmbedingDict[nbrAtom.GetIdx()])) :
447                            candidates.append(FragmentExtensionByBond(atomInd, bond.GetBondOrder(), reverseEmbedingDict[nbrAtom.GetIdx()], [embeding]))
448                    else:
449                        candidates.append(FragmentExtensionByBondAtom(atomInd, bond.GetBondOrder(), [nbrAtom], [embeding]))
450        #candidates.sort(lambda a,b:cmp(a[0],b[0]) or cmp(a[1].GetAtomicNum(),b[1].GetAtomicNum()))
451        return candidates
452
453    def GetRingCandidatesFromRing(self, ring, embeding):
454        candidates=[]
455        reverseEmbedingDict=embeding.GetReverseDict()
456        for candidateRing in self.FilterRings(self.miner.rings[embeding.molecule], embeding):
457            for atom in ring.ringAtoms:
458                if candidateRing.IsMember(embeding.GetEmbededAtom(atom.atomIndex)):
459                    candidates.append(FragmentExtensionByRing([candidateRing], [embeding]))
460                    break
461        return candidates                             
462   
463    def GetCandidatesFromEmbeding(self, embeding):
464        candidates=[]
465        if len(self.atoms)==1 and self.miner.addWholeRings:
466            atom=embeding.GetEmbededAtom(self.atoms[0].GetIdx())
467            if atom.IsInRing():
468                for ring in embeding.molecule.rings:
469                    if ring.IsMember(atom):
470                        candidates.append(FragmentExtensionByRing([ring],[embeding]))
471                return candidates
472           
473        if debug: print "Parsing candidate bonds"
474        for atom in self.GetCandidateAtoms():
475            candidates.extend(self.GetCandidatesFromAtom(atom.atomIndex, embeding))
476       
477        if debug: print "Parsing candidate rings"
478        for ring in self.GetCandidateRings():
479            candidates.extend(self.GetRingCandidatesFromRing(ring, embeding))           
480        return candidates       
481       
482    def Extend(self):
483        if debug :print "Extending"
484        candidates=[]
485        for embeding in self.embedings:
486            c=self.GetCandidatesFromEmbeding(embeding)
487            candidates.extend(c)
488            #embedingDict[embeding]=c
489           
490        #group equivalent candidates
491        if debug: print "Grouping candidates"
492        groups=[]
493        for extension in candidates:
494            for ext in groups:
495                if ext.__class__==extension.__class__ and ext.IsEquivalent(extension):
496                    ext.MergeFrom(extension)
497                    break
498            else:
499                groups.append(extension)
500 
501        #generate new fragments
502        if debug: print "Generating new fragments"
503        newFragments=[]
504        for extension in groups:
505            #set_trace()
506            f=deepcopy(self)
507            extension.Extend(f)
508            newFragments.append(f)
509        return newFragments
510   
511    def ToOBMol(self):
512        atomCache={}
513        mol=OBMol()
514        mol.BeginModify()
515        for sourceAtom in self.atoms:
516            atom=mol.NewAtom()
517            atom.SetAtomicNum(sourceAtom.GetAtomicNum())
518            if sourceAtom.IsAromatic():
519                atom.SetAromatic()
520##                atom.SetSpinMultiplicity(2)
521            atomCache[sourceAtom]=atom
522        for sourceBond in self.bonds:
523            mol.AddBond(atomCache[sourceBond.atom1].GetIdx(), atomCache[sourceBond.atom2].GetIdx(), sourceBond.GetBondOrder())
524##        mol.SetAromaticPerceived()
525        mol.AssignSpinMultiplicity()
526##        mol.UnsetAromaticPerceived()
527        mol.EndModify()
528        return mol
529
530    def ToSmiles(self):
531        writer=OBConversion()
532        writer.SetInAndOutFormats("smi", "smi")
533        return writer.WriteString(self.ToOBMol()).strip()
534   
535    def ToCannonicalSmiles(self):
536        atomCache={}
537        mol=OBMol()
538        for sourceAtom in self.atoms:
539            atom=mol.NewAtom()
540            atom.SetAtomicNum(sourceAtom.GetAtomicNum())
541            if sourceAtom.IsAromatic():
542                atom.SetAromatic()
543                atom.SetSpinMultiplicity(2)
544            atomCache[sourceAtom]=atom
545        for sourceBond in self.bonds:
546            mol.AddBond(atomCache[sourceBond.atom1].GetIdx(), atomCache[sourceBond.atom2].GetIdx(), sourceBond.GetBondOrder())
547        writer=OBConversion()
548        writer.SetInAndOutFormats("smi", "can")
549        return writer.WriteString(mol).strip()
550
551    def Support(self, activeSet=None):
552        activeSet=self.miner.activeSet if activeSet==None else activeSet
553        uniqueMolecules=set()
554        for embeding in self.embedings:
555            if embeding.molecule in activeSet:
556                uniqueMolecules.add(embeding.molecule)
557##        s=set(filter(lambda mol:self.ContainedIn(mol), self.miner.GetAllMolecules()))
558##        if len(s) != len(uniqueMolecules):
559##            writer=OBConversion()
560##            writer.SetInAndOutFormats("smi", "smi")
561##            print "\n",self.ToSmiles()
562##            for m in uniqueMolecules: print writer.WriteString(m).strip()
563##            for m in s: print writer.WriteString(m).strip()
564       
565        return float(len(uniqueMolecules))/float(len(activeSet) or 1)
566
567    def OcurrencesIn(self, molecule):   
568        pattern=OBSmartsPattern()
569        pattern.Init(self.ToSmiles())
570        return pattern.Match(molecule)
571   
572    def ContainedIn(self, molecule):
573        return bool(self.OcurrencesIn(molecule))
574           
575class FragmentMiner(object):
576    """A class for finding frequent molecular fragments
577    Attributes:
578        active      : list of smiles codes of active molecules
579        inactive    : list of smiles codes of inactive molecules
580        minSupport  : minimum frequency in the active set of the fragments to search for
581        maxSupport  : maximum frequency in the inactive set of the fragments to search for
582        addWholeRings : if True rings will be added as a whole rather then atom by atom
583        canonicalPruning : if True a cache of all cannonical codes of all fragments will be kept to avoid
584                    redundant search
585        findClosed  : finds only fragments that are not sub-structures of any other fragment with the same support (default: True)
586    Example:
587    >>> miner = FragmentMiner(active = ["CC(C=N)=O", "c1ccccc1C=O", "SCC(N)O"], inactive = [], minSupport = 0.6)
588    >>> for fragment in miner.Search():
589    ...     print fragment.ToSmiles() , "Support: %.3f" %fragment.Support()
590    """
591    def __init__(self, active, inactive=[], minSupport=0.2, maxSupport=0.2, addWholeRings=True, canonicalPruning=True, findClosed=True):
592        self.active=filter(lambda m:m, map(self.LoadMolecules, active))
593        self.inactive=filter(lambda m:m, map(self.LoadMolecules, inactive))
594        self.minSupport=minSupport
595        self.maxSupport=maxSupport
596        self.rings={}
597        self.atomCount={}
598        self.findClosed=findClosed
599        self.addWholeRings=addWholeRings
600        self.canonicalPruning=canonicalPruning
601        self.canonicalPruningSet={}
[953]602        self.loader=OBConversion()
603        self.loader.SetInAndOutFormats("smi","smi")
[597]604
605    def LoadMolecules(self, smiles):
606        mol=LoadMolFromSmiles(smiles)
[644]607##        if mol:
608##            mol.StripSalts()
[597]609        return mol
610       
611    def GetAllMolecules(self):
612        return self.active+self.inactive     
613
614    def Initialize(self):
615        """Initializes the search"""
616        self.initialFragments=[]
617        self.rings={}
618        self.atomCount={}
619        self.canonicalPruningSet={}
620        candidates=[]
621        ringCandidates=[]
622        for mol in self.GetAllMolecules():
623            mol.rings=self.rings[mol]=list(mol.GetSSSR())
624            for ring in mol.rings:
625                ring.fingerprint=set([mol.GetAtom(i).GetAtomicNum() for i in ring._path])
626##            for ring in mol.rings:
627##                candidates.append(FragmentExtensionByRing([ring],[Embeding(molecule=mol)]))
628       
629        for mol in self.GetAllMolecules():
630            for atom in OBMolAtomIter(mol):
631                self.atomCount[atom.GetAtomicNum()]= self.atomCount[atom.GetAtomicNum()]+1 if  atom.GetAtomicNum() in self.atomCount else 1
632                if not (atom.GetAtomicNum()==6 and atom.IsInRing()):
633                    candidates.append(FragmentExtensionByAtom(atom.GetAtomicNum(), atom.IsAromatic(), [atom.GetIdx()], [Embeding(molecule=mol)]))
634        groups=[]
635        candidates.sort(lambda a,b: cmp(self.atomCount[a.atomicNum], self.atomCount[b.atomicNum]))
636        for extension in candidates:
637            for ext in groups:
638                if type(ext)==type(extension) and ext.IsEquivalent(extension):
639                    ext.MergeFrom(extension)
640                    break
641            else:
642                groups.append(extension)
643        lst=self.atomCount.items()
644        lst.sort(lambda a,b:cmp(a[1], b[1]))
645        self.initialFragments=[]
646        lst=[t[0] for t in lst]
647        for extension in groups:
648            if type(extension)==FragmentExtensionByAtom:
649                f=Fragment(miner=self, excludeAtomList=lst[:lst.index(extension.atomicNum)])
650            else:
651                f=Fragment(miner=self)
652            extension.Extend(f)
653            self.initialFragments.append(f)
654##        self.initialFragments.reverse()
655##        excludeList=[]
656##        for atom, c in lst:
657##            f=Fragment(miner=self, excludeAtomList=excludeList)
658##            f.InitializeFragment(atom)
659##            self.initialFragments.append(f)
660##            excludeList.append(atom)
661        self.activeSet=set(self.active)
662        self.inactiveSet=set(self.inactive)
663           
664    def TraverseTree(self, fragment):
665        if self.canonicalPruning:
666            codeWord=fragment.ToCannonicalSmiles()
667            if codeWord in self.canonicalPruningSet:
668                return self.canonicalPruningSet[codeWord]
669            else:
670                self.canonicalPruningSet[codeWord]=fragment.Support(self.activeSet)
671        extended=fragment.Extend()
672        extended=filter(lambda f:f.Support(self.activeSet)>=self.minSupport, extended)
673        superStructSupport=[]
674        for frag in extended:
675            #print self.loader.WriteString(frag.ToOBMol())
676            superStructSupport.append(self.TraverseTree(frag))
677        support=fragment.Support(self.activeSet)
678        if support>=self.minSupport and fragment.Support(self.inactiveSet)<=self.maxSupport:
679            if not self.findClosed or (support not in superStructSupport):
680                print fragment.ToSmiles().strip()+" %.2f %.2f" % (support, fragment.Support(self.inactiveSet))
681                self.foundFragments.append(fragment)
682        return support
683
684    def Search(self):
685        """Runs the search and returns the found fragments"""
686        self.Initialize()
687##        set_trace()
688        self.foundFragments=[]
689        for fragment in self.initialFragments:
690            self.TraverseTree(fragment)
691        #self.foundFragments=filter(lambda f:f.Support(self.inactive)<=self.maxSupport, self.foundFragments)
692        return self.foundFragments
693
694    def TraverseTreeIterator(self, fragment):
695        if self.canonicalPruning:
696            codeWord=fragment.ToCannonicalSmiles()
697            if codeWord in self.canonicalPruningSet:
698                raise StopIteration
699            else:
700                self.canonicalPruningSet[codeWord]=fragment.Support(self.activeSet)
701        extended=fragment.Extend()
702        extended=filter(lambda f:f.Support(self.activeSet)>=self.minSupport, extended)
703        superStructSupport=[]
704        for frag in extended:
705            #print self.loader.WriteString(frag.ToOBMol())
706            iter=self.TraverseTreeIterator(frag)
707            try:
708                while True:
709                    f=iter.next()
710                    superStructSupport.append(f.Support(self.active))
711                    yield f
712            except StopIteration:
713                pass
714               
715            superStructSupport.append(self.TraverseTree(frag))
716        support=fragment.Support(self.activeSet)
717        if support>=self.minSupport and fragment.Support(self.inactiveSet)<=self.maxSupport:
718            if not self.findClosed or (support not in superStructSupport):
719                #print fragment.ToSmiles().strip()+" %.2f %.2f" % (support, fragment.Support(self.inactiveSet))
720                #self.foundFragments.append(fragment)
721                yield fragment
722        #return support
723
724    def SearchIterator(self):
725        """Runs the search and returns the found fragments one by one"""
726        self.Initialize()
727##        set_trace()
728        self.foundFragments=[]
729        for fragment in self.initialFragments:
730            iter=self.TraverseTreeIterator(fragment)
731            try:
732                while True:
733                    yield iter.next()
734            except StopIteration:
735                pass
736        #self.foundFragments=filter(lambda f:f.Support(self.inactive)<=self.maxSupport, self.foundFragments)
737        #return self.foundFragments   
738   
739def LoadMolFromSmiles(smiles):
740    """Returns an OBMol construcetd from an SMILES code"""
[644]741    smiles = sorted(smiles.split("."), key=len)[-1] ## Strip salts
[597]742    mol=OBMol()
743    loader=OBConversion()
744    loader.SetInAndOutFormats("smi","smi")
745    if not loader.ReadString(mol, smiles):
746        return None
747    mol.smilesCode=smiles
748    return mol
749   
750class Fragmenter(object):
751    """An object that is used to fragment an ExampleTable
752    Attributes:
753        minSupport  : minimum frequency in the active set of the fragments to search for (default: 0.2)
754        maxSupport  : maximum frequency in the inactive set of the fragments to search for (default: 0.2)
755        findClosed  : finds only fragments that are not sub-structures of any other fragment with the same support (default: True)
756    Example:
757    >>> fragmenter=Fragmenter(minSupport=0.1, maxSupport=0.05)
758    >>> data, fragments=fragmenter(data, "SMILES", lambda ex:ex.getclass())
759    """
760    def __init__(self, minSupport=0.2, maxSupport=0.2, canonicalPruning=True, findClosed=True):
761        self.minSupport=minSupport
762        self.maxSupport=maxSupport
763        self.canonicalPruning=canonicalPruning
764        self.findClosed=findClosed
[887]765    def __call__(self, data, smilesAttr=None, activeFunc=lambda e:True, useCannonicalFragments=True):
[597]766        """Takes a data-set, and runs the FragmentMiner on it. Returns a new data-set and the fragments.
767        The new data-set contains new attributes that represent the presence of a fragment that was found.
768        Arguments:
769            data        : the dataset
770            smilesAttr  : the attribute in the data that contains the SMILES codes
771            activeFunc  : a function that takes an example from the data-set and returns True if the example should be
772                    considered as active (if none is provided all examples are considered active)
[887]773            useCannonicalFragments  : if True use cannonical smiles for fragment attribute names else use regular smiles (defaut: True)
[597]774        """
775        if not smilesAttr:
776            smilesAttr=self.FindSmilesAttr(data)
777        active=filter(lambda s:s, [str(e[smilesAttr]) for e in data if activeFunc(e)])
778        inactive=filter(lambda s:s, [str(e[smilesAttr]) for e in data if not activeFunc(e)])
779       
780        miner=FragmentMiner(active, inactive, self.minSupport, self.maxSupport, canonicalPruning=self.canonicalPruning, findClosed=self.findClosed)
781        self.fragments=fragments=miner.Search()
[887]782        fragVars=[orange.EnumVariable(frag.ToCannonicalSmiles() if useCannonicalFragments else frag.ToSmiles(), values=["0", "1"]) for frag in fragments]
[597]783        smilesInFragments=dict([(fragment, set([embeding.molecule.smilesCode for embeding in fragment.embedings]) ) for fragment in fragments])
784        from functools import partial
785        def getVal(var, fragment, smilesAttr, example, returnWhat):
786            mol=LoadMolFromSmiles(str(example[smilesAttr]))
787##            print "GetVal"
[885]788            return (fragment.ContainedIn(mol) and var(1) or var(0)) if mol else var(None)
[597]789        for var, frag in zip(fragVars, fragments):
790            var.getValueFrom=partial(getVal,var, frag, smilesAttr)
791        vars=data.domain.attributes+fragVars+(data.domain.classVar and [data.domain.classVar] or [])
792        domain=orange.Domain(vars, data.domain.classVar and 1 or 0)
793        domain.addmetas(data.domain.getmetas())
794        table=orange.ExampleTable(domain)
795        for e in data:
[886]796            vals=[e[attr] for attr in data.domain.attributes]+[1 if str(e[smilesAttr]) in smilesInFragments[fragment] else (None if e[smilesAttr].isSpecial() else 0) for fragment in fragments]
[597]797            vals=vals + [e.getclass()] if data.domain.classVar else vals
798            ex=orange.Example(domain, vals)
799            for key, val in e.getmetas().items():
800                ex[key]=val
801            table.append(ex)
802        return table, fragments
803    def FindSmilesAttr(self, data):
804        data=data.select(orange.MakeRandomIndices2(data, min(20, len(data))))
805        stringVars=filter(lambda var:type(var)==orange.StringVariable, data.domain.attributes+data.domain.getmetas().values())
806        count=dict.fromkeys(stringVars, 0)
807        for example in data:
808            for var in stringVars:
809                if LoadMolFromSmiles(str(example[var])):
810                    count[var]+=1
811        count=count.items()
812        count.sort(lambda a,b:cmp(a[1], b[1]))
813        return count[-1][0]
814           
[1632]815from Orange.orng import orngSVM
[597]816class FragmentBasedLearner(orange.Learner):
817    """A learner wrapper class that first runs the molecular fragmentation on the data.
818    Attributes:
819        smilesAttr  : Attribute in the data that contains the smiles codes (if none is provided it will try to make a smart guess)
820        learner     : learner that will be used to actualy learn on the fragmented data (default: orngSVM.SVMLearner)
821        minSupport  : minimum frequency in the active set of the fragments to search for
822        maxSupport  : maximum frequency in the inactive set of the fragments to search for
823        activeFunc  : a function that takes an example from the learning data-set and returns True if the example should be
824                    considered as active (if none is provided all examples are considered active)
825        findClosed  : finds only fragments that are not sub-structures of any other fragment with the same support (default: True)
826    """
827    def __new__(cls, data=None, weights=0, **kwds):
828        learner=orange.Learner.__new__(cls, **kwds)
829        if data:
830            learner.__init__(**kwds)
831            return learner(data)
832        else:
833            return learner
834    def __init__(self, learner=orngSVM.SVMLearner(probability=True), name="FragmentBasedLearner",
835                 minSupport=0.2, maxSupport=0.2, smilesAttr=None, findClosed=True, activeFunc=lambda e:True):
836        self.name=name
837        self.learner=learner
838        self.minSupport=minSupport
839        self.smilesAttr=smilesAttr
840        self.activeFunc=activeFunc
841        self.maxSupport=maxSupport
842        self.findClosed=findClosed
843    def __call__(self, data, weight=0):
844        fragmenter=Fragmenter(minSupport=self.minSupport, maxSupport=self.maxSupport, findClosed=self.findClosed)
845        data, fragments=fragmenter(data, self.smilesAttr, self.activeFunc)
846        return FragmentBasedClassifier(self.learner(data), data.domain)
847
848class FragmentBasedClassifier(object):
849    def __init__(self, classifier, domain):
850        self.classifier=classifier
851        self.domain=domain
852    def __call__(self, example, getBoth=orange.GetValue):
853        example=orange.Example(self.domain, example)
854        return self.classifier(example, getBoth)
855
856def Count(smiles, fragment):
857    mols=filter(lambda m:m, map(LoadMolFromSmiles, smiles))
[644]858##    for mol in mols: mol.StripSalts()
[597]859    pattern=OBSmartsPattern()
860    pattern.Init(fragment)
861    return len(filter(lambda m:pattern.Match(m, True), mols))
862
863def ContaindIn(smiles, fragment):
864    mol=LoadMolFromSmiles(smiles)
865    pattern=OBSmartsPattern()
866    pattern.Init(fragment)
867    return bool(pattern.Match(mol))
[598]868
[1320]869
870"""
871Viasualizations
872===============
873
874"""
[598]875try:
[606]876    from oasa.svg_out import svg_out as _svg_out
[1320]877    from oasa.pybel_bridge import PybelConverter
878    from oasa.coords_generator import coords_generator
879    from oasa import dom_extensions
[598]880except ImportError:
[1320]881    try:
882        from bkchem.oasa.oasa.svg_out import svg_out as _svg_out
883        from bkchem.oasa.oasa.pybel_bridge import PybelConverter
884        from bkchem.oasa.oasa.coords_generator import coords_generator
885        from bkchem.oasa.oasa import dom_extensions
886    except ImportError:
887        class _svg_out(object):
888            def __init__(self, *args, **kwargs):
889                raise ImportError("OASA library not found")
[695]890
891try:
892    from oasa.cairo_out import cairo_out as _cairo_out
893except ImportError:
[1320]894    try:
895        from bkchem.oasa.oasa.cairo_out import cairo_out as _cairo_out
896    except ImportError:
897        class _cairo_out(object):
898            def __init__(self, *args, **kwargs):
899                raise ImportError("OASA library not found")
[598]900
[601]901class cairo_out(_cairo_out):
[598]902    atom_colors = {"A":(0, 0, 0), "R":(255, 0, 0)}
903    def __init__(self, **kwargs):
[601]904        _cairo_out.__init__(self, **kwargs)
[598]905        self.atom_colors.update({"A":(0, 0, 0), "R":(255, 0, 0)})
906
907    def _draw_edge(self, e):
908        v1, v2 = e.vertices
909        s1, s2 = v1.symbol, v2.symbol
910        matched = v1 in self.ob_matched and v2 in self.ob_matched
[601]911        ## 'A' and 'R' are used because they are in oasa periodic table and indicate any atom or any atom except hydrogen
[598]912        v1.symbol = "R" if matched else "A"
913        v2.symbol = "R" if matched else "A"
[599]914        if matched:
[601]915            ## Set the default color argument for _draw_line to red
[599]916            self._draw_line.im_func.func_defaults = self._draw_line.im_func.func_defaults[:2]+ ((255, 0, 0),)
[601]917        _cairo_out._draw_edge(self, e)
[599]918        if matched:
[601]919            ## Restore the default color
[599]920            self._draw_line.im_func.func_defaults = self._draw_line.im_func.func_defaults[:2]+ ((0, 0, 0),)
[598]921        v1.symbol, v2.symbol = s1, s2
922
923    def _draw_vertex(self, v):
924        if v in self.ob_matched:
925            color = self.atom_colors.get(v.symbol, None)
926            self.atom_colors[v.symbol] = (255, 0, 0)
[601]927            _cairo_out._draw_vertex(self, v)
928            ## Restore the atom colors
[598]929            if color:
930                self.atom_colors[v.symbol] = color
931            else:
932                del self.atom_colors[v.symbol]
933        else:
[601]934            _cairo_out._draw_vertex(self, v)
[598]935
936    def mol_to_cairo(self, molSmiles, fragSmiles=None, file="mol.png"):
[601]937        import pybel
[606]938        molSmiles = max(molSmiles.split("."), key = lambda s:len(s))
939        mol = pybel.readstring("smi", molSmiles)
[601]940       
[598]941        if fragSmiles:
[601]942            pattern=pybel.Smarts(fragSmiles)
943            matches = pattern.findall(mol)
[598]944            self.ob_matched = reduce(set.union, matches if matches != False else [], set())
945        else:
946            self.ob_matched = set()
[1320]947#        from oasa.pybel_bridge import PybelConverter
[601]948        o_mol, idx2oa = PybelConverter.pybel_to_oasa_molecule_with_atom_map(mol)
[1320]949#        from oasa.coords_generator import coords_generator
[598]950        gen = coords_generator(30)
951        gen.calculate_coords(o_mol, bond_length=30, force=True)
952##        self.oasa2obidx = dict([(value, key) for key, value in ix2oa.items()])
953        self.ob_matched = [idx2oa[i] for i in self.ob_matched]
954        self.color_bonds = True
[606]955        return _cairo_out.mol_to_cairo(self, o_mol, file)
[601]956
[606]957import xml.dom.minidom as dom
958class svg_out(_svg_out):
959    def __init__(self, **kwargs):
960        _svg_out.__init__(self, **kwargs)
961        self._curr_color = "rgb(0%, 0%, 0%)"
962
963    def _draw_edge(self, e):
964        v1, v2 = e.vertices
965        matched = v1 in self.ob_matched and v2 in self.ob_matched
966        if matched:
967            self._curr_color = "rgb(255, 0, 0)"
968        else:
969            self._curr_color = "rgb(0, 0, 0)"
970
971        _svg_out._draw_edge(self, e)
972
973    def _draw_vertex(self, v):
974        if v in self.ob_matched:
975            self._curr_color = "rgb(255, 0, 0)"
976        else:
977            self._curr_color = "rgb(0, 0, 0)"
978        _svg_out._draw_vertex(self, v)
979
980
981    def _draw_line( self, parent, start, end, line_width=1, capstyle=""):
982        x1, y1 = start
983        x2, y2 = end
[1320]984#        from oasa import dom_extensions
[606]985        line = dom_extensions.elementUnder( parent, 'line',
986                                            (( 'x1', str( x1)),
987                                             ( 'y1', str( y1)),
988                                             ( 'x2', str( x2)),
989                                             ( 'y2', str( y2)),
990                                             ( 'stroke', self._curr_color)))
991
992
993    def _draw_text( self, parent, xy, text, font_name="Arial", font_size=16):
994        x, y = xy
[1320]995#        from oasa import dom_extensions
[606]996        dom_extensions.textOnlyElementUnder( parent, "text", text,
997                                             (( "x", str( x)),
998                                              ( "y", str( y)),
999                                              ( "font-family", font_name),
1000                                              ( "font-size", str( font_size)),
1001                                              ( 'stroke', self._curr_color)))
1002
1003    def mol_to_svg(self, molSmiles, fragSmiles=None):
1004        import pybel
1005        molSmiles = max(molSmiles.split("."), key = lambda s:len(s))
1006        mol = pybel.readstring("smi", molSmiles)
1007       
1008        if fragSmiles:
1009            pattern=pybel.Smarts(fragSmiles)
1010            matches = pattern.findall(mol)
1011            self.ob_matched = reduce(set.union, matches if matches != False else [], set())
1012        else:
1013            self.ob_matched = set()
[1320]1014#        from oasa.pybel_bridge import PybelConverter
[606]1015        o_mol, idx2oa = PybelConverter.pybel_to_oasa_molecule_with_atom_map(mol)
[1320]1016#        from oasa.coords_generator import coords_generator
[606]1017        gen = coords_generator(30)
1018        gen.calculate_coords(o_mol, bond_length=30, force=True)
1019##        self.oasa2obidx = dict([(value, key) for key, value in ix2oa.items()])
1020        self.ob_matched = [idx2oa[i] for i in self.ob_matched]
1021        return _svg_out.mol_to_svg(self, o_mol)   
1022       
[601]1023def mol_to_png(molSmiles, fragSmiles=None, file="mol.png", size=None):
1024    d = cairo_out()
[606]1025    d.mol_to_cairo(molSmiles, fragSmiles)
1026    if size != None:
1027        import Image
1028        image = Image.open(file)
1029        width, height = image.size
1030        resize_factor = float(size) / max(width, height)
1031        width, height = int(width * resize_factor), int(height * resize_factor)
1032        image = image.resize((width, height), Image.ANTIALIAS)
1033        newimage = Image.new("RGB", size=(size, size), color=(255, 255, 255))
1034        newimage.paste(image, (size - width, size - height))
1035        newimage.save(file)
1036
1037def mol_to_svg(molSmiles, fragSmiles=None, file="mol.svg"):
1038    c = svg_out()
1039    tree = c.mol_to_svg(molSmiles, fragSmiles)
1040    if type(file) == str:
1041        file = open(file, "w")
1042    file.write( tree.toprettyxml())
[601]1043   
[598]1044def _test_draw():
[601]1045##    mol_to_png("CC1CC2C(CC(=O)O2)OC3CC4C(CC(C5C(O4)CC=CCC6C(O5)CC=CC7C(O6)CCCC8C(O7)(CC9C(O8)CC2C(O9)C(CC(O2)CC(=C)C=O)O)C)C)OC3(C1)C", "Scc", "mol1.png")
[606]1046##    mol_to_png("CN(C)CCCN1c2ccccc2Sc3c1cc(cc3)C(F)(F)F", "Scc", "mol1.png")
1047##    mol_to_png("CN(C)CCCN1c2ccccc2Sc3c1cc(cc3)C(F)(F)F", "NCC", "mol2.png")
1048##    mol_to_png("CN(C)CCCN1c2ccccc2Sc3c1cc(cc3)C(F)(F)F", "Ncc", "mol3.png")
1049##    mol_to_png("CN(C)CCCN1c2ccccc2Sc3c1cc(cc3)C(F)(F)F", "NCCC", "mol3.png")
1050    mol_to_svg("CN(C)CCCN1c2ccccc2Sc3c1cc(cc3)C(F)(F)F", "Scc", "mol1.svg")
1051    mol_to_svg("CN(C)CCCN1c2ccccc2Sc3c1cc(cc3)C(F)(F)F", "Ncc", "mol2.svg")
1052    mol_to_svg("CN(C)CCCN1c2ccccc2Sc3c1cc(cc3)C(F)(F)F", "NCC", "mol3.svg")
[597]1053   
[598]1054def _test():
[597]1055    import orange
1056    d=orange.ExampleTable("E:\chem\mutagen_raw.tab")
1057    active=[str(e["SMILES"]) for e in d if str(e[-1])=="1"]
1058    inactive=[str(e["SMILES"]) for e in d if str(e[-1])=="0"]
1059##    d=orange.ExampleTable("E:\PCLedit_s.tab")
1060##    active=[str(e["SMILES"]) for e in d if not e["SMILES"].isSpecial()][:100]
1061##    print active
1062##    inactive=[]
1063##    active=["NC(C)C(=O)O", "NC(CS)C(=O)O", "NC(CO)C(=O)O"]
1064##    active=["CCS(O)(O)N", "CCS(O)(C)N", "CS(O)(C)N", "CCS(=N)N", "CS(=N)N", "CS(=N)O"]
1065##    active=["NC(S)c1ccccc1","NCC1=CC=CC=C1", "NCC1C=CC=CC=1", "c1ccccc1C(N)C(=S(O)C)c2ccccc2"]
1066##    active=["c1ccccc1C(N)C(=S(O)C)c2ccccc2"]
1067##    active=["O=C1C=CC(=O)C=C1","O=C1CCCCCN1"]
1068##    active=["C1SC2CCN2C1C(=O)"]
1069##    active=["CCCCCCc1ccc(O)cc1O","Nc1ccc(O)c(N)c1", "Cc1cc(C)c(N)cc1C", "CN(C)C(=S)S[Zn]SC(=S)N(C)C", "NC(=O)N(CCO)N=O"]
1070##    active=["CC(C)CCCC(C)C1CCC2C3CC=C4CC(CCC4(C)C3CCC12C)OC(=O)Cc5ccc(cc5)N(CCCl)CCCl","Cc1cc(C)c(N=Nc2c(O)c(cc3cc(ccc23)S(=O)(=O)O)S(=O)(=O)O)cc1CCNNCc1ccc(cc1)C(=O)NC(C)C","CC(C)(Oc1ccc(cc1)C2CCCc3ccccc23)C(=O)O","CCn1cc(C(=O)O)c(=O)c2ccc(C)nc12"]
1071    active=["CN(C)CCCN1c2ccccc2Sc3c1cc(cc3)C(F)(F)F", "CN(C)CCCN1c2ccccc2Sc3c1cc(cc3)Cl","CN1CCCCC1CCN2c3ccccc3Sc4c2cc(cc4)SC","c1ccc2c(c1)N(c3cc(ccc3S2)Cl)CCCN4CCN(CC4)CCO",
1072            "CN1CCN(CC1)CCCN2c3ccccc3Sc4c2cc(cc4)S(=O)(=O)N(C)C", "CN1CCC(=C2c3ccccc3Sc4c2cccc4)CC1", "[U]-C-S-P"]
1073##    active=["CN(C)CCCN1c2ccccc2Sc3c1cc(cc3)C(F)(F)F", "CN(C)CCCN1c2ccccc2Sc3c1cc(cc3)Cl","CN1CCCCC1CCN2c3ccccc3Sc4c2cc(cc4)SC"]
1074##    active=["CN(C)CCCN1c2ccccc2Sc3c1cc(cc3)C(F)(F)F", "CN(C)CCCN1c2ccccc2Sc3c1cc(cc3)Cl"]
1075##    set_trace()
1076    miner=FragmentMiner(active, inactive[:0], minSupport=0.1, maxSupport=0.1, addWholeRings=True, canonicalPruning=True)
1077    fragments=miner.Search()
1078##    for f in fragments:
1079##        print f.ToSmiles()
1080
[598]1081def _test1():
[597]1082    import orange
1083    data=orange.ExampleTable("E:\chem\mutagen_raw.tab")
1084##    data=orange.ExampleTable("E:\chem\smiles.tab")
1085    fragmenter=Fragmenter(minSupport=0.02, maxSupport=0.1, canonicalPruning=True)
1086##    set_trace()
1087    data, fragments1=fragmenter(data, "SMILES") #, lambda e:str(e[-1])=="1")
1088##    data, fragments2=fragmenter(data, "SMILES", lambda e:str(e[-1])=="0")
1089    data.save("E:\chem\mutagen_raw_frag.tab")
1090   
1091if __name__=="__main__":
1092    import time
1093    sTime=time.clock()
[598]1094    _test_draw()
[597]1095    print time.clock()-sTime
1096
Note: See TracBrowser for help on using the repository browser.