source: orange-bioinformatics/_bioinformatics/obiMeSH.py @ 1637:dde4ef52193b

Revision 1637:dde4ef52193b, 30.9 KB checked in by mitar, 2 years ago (diff)

Removing __file__ usage.

Line 
1"""
2Module for browsing any analyzing sets annotated with MeSH ontology.
3"""
4
5from xml.sax import make_parser
6from xml.dom import minidom
7from xml.sax.handler import ContentHandler
8from math import log,exp
9from urllib import urlopen
10from sgmllib import SGMLParser
11import os.path
12
13import orange
14from Orange.orng import orngServerFiles
15
16FUZZYMETAID = -12
17
18class obiMeSH(object):
19    def __init__(self):
20        """
21        Function will initialize obiMeSH module and (in case of no MeSH ontology data) download necessary data using orngServerFiles.
22        """
23        #print "SVN version"
24        # we check if all files from MeSH directory are localy present. if not, download them
25        d = orngServerFiles.ServerFiles()
26        f = d.listfiles('MeSH')
27        l = orngServerFiles.listfiles('MeSH')
28        for i in f:
29            if not i in l:
30                orngServerFiles.download('MeSH',i)
31                print 'Downloading MeSH ontology and chemical annotation. Please wait.', i
32       
33        self.reference = None
34        self.cluster = None
35        self.ratio = 1
36        self.statistics = None
37        self.calculated = False
38       
39        self.ref_att = "Unknown"
40        self.clu_att = "Unknown"
41        self.solo_att = "Unknown"
42       
43        #we calculate log(i!)
44        self.lookup = [0]
45        for i in range(1, 8000):
46            self.lookup.append(self.lookup[-1] + log(i))
47        self.dataLoaded = self.__loadOntologyFromDisk()
48
49    def expandToFuzzyExamples(self, examples, att, a, b):
50        """
51        Function will return new 'fuzzy' example table. Every example from the input table will get two additional meta attributes ('fuzzy set' and 'u') \
52        based on 'a' and 'b' threshold (lower and higher) and attribute 'att'. Attribute 'fuzzy set' indicates name of the fuzzy set while atribute 'u' \
53        reflects example's degree of membership to particular fuzzy set. Note that input examples with values of 'att' lying on the (a,b) will be expanded \
54        into two fuzzy examples.
55        """
56        mu = orange.FloatVariable("u")
57        mv = orange.StringVariable("fuzzy set")
58        examples.domain.addmeta(FUZZYMETAID, mu)
59        examples.domain.addmeta(FUZZYMETAID-1, mv)
60        newexamples = []
61        for j in range(0,len(examples)):
62            i = examples[j]
63            v = float(i[att])
64            if v > a and v < b:  # we have to expand this example
65                newexamples.append(i)
66                i["fuzzy set"] = 'yes'
67                i["u"] = (v-a)/(b-a)
68                examples.append(i)
69                examples[-1]["fuzzy set"] = "no"
70                examples[-1]["u"] = (b-v)/(b-a)
71            else:
72                if v > a:  #        u(yes) = 1.0
73                    i["fuzzy set"] = 'yes'
74                    i["u"] = 1.0
75                else: #     u(no) = 1.0
76                    i["fuzzy set"] = 'no'
77                    i["u"] = 1.0
78        return examples
79
80    def findSubset(self, examples, meshTerms, callback=None, MeSHtype='term'):
81        """
82        Function findSubset will return new example table containing examples which are annotated with at least one MeSH term from list 'meshTerms'.
83        """
84        # clone
85        newdata = orange.ExampleTable(examples.domain)
86        self.solo_att = self.__findMeshAttribute(examples)
87        ids = list()
88        l = len(examples)
89        c = 0.0
90        meshTerms = list(set(meshTerms))
91        # we couldn't find any mesh attribute
92        if self.solo_att == "Unknown":
93            return newdata
94        if MeSHtype == 'term':
95            for i in meshTerms:
96                ids.extend(self.toID[i])
97        else:
98            ids = meshTerms
99        for e in examples:
100            try:
101                if callback and c%10 == 0:
102                    callback(int(c*100.0/l))
103                c = c + 1.0
104                ends = eval(e[self.solo_att].value)
105            except SyntaxError:
106                #print "Error in parsing ", e[self.solo_att].value
107                continue
108            endids = list()
109            for i in ends:
110                if self.toID.has_key(i):
111                    endids.extend(self.toID[i])
112            #allnodes = self.__findParents(endids)
113            allnodes = endids
114            # calculate intersection
115            isOk = False
116            for i in allnodes:
117                if ids.count(i) > 0:
118                    isOk = True
119                    break
120
121            if isOk:      # intersection between example mesh terms and observed term group is None
122                newdata.append(e)
123        return newdata
124
125    def findTerms(self, ids, idType="cid", callback=None):
126        """
127        Function findTerms returns a dictionary containing annotated items from the list 'ids'.
128        """
129        ret = dict()
130        if(not self.dataLoaded):
131            print "Annotation and ontology has never been loaded!"
132            return ret
133
134        if idType=="cid": 
135            for i in ids:
136                if self.fromCID.has_key(int(i)): # maybe it is on our database
137                    ret[int(i)] = self.fromCID[int(i)]
138                else:   # no it is not, let's go on the internet - use pubChemAPI
139                    r = pubChemAPI()
140                    l = r.getMeSHterms(i)
141                    ca = r.getPharmActionList(i)
142                    l.extend(ca) # we extend it with Pharmacological actions
143                    #l = self.findTermsForCID(i)
144                    self.fromCID[int(i)] = l
145                    ret[int(i)] = l
146                   
147                    #if len(l)>0:
148                    #   # if we found something lets save it to a file
149                    #   TODO: Never use __file__ to access data, but pkg_resources so that it works also in eggs
150                    #   __dataPath = os.path.join(os.path.dirname(__file__), self.path)
151                    #   fileHandle = open(os.path.join(__dataPath,'cid-annotation.dat'), 'a')
152                    #   for s in l:
153                    #       fileHandle.write('\n' + str(i) + ';' + s )
154                    #   fileHandle.close()
155            return ret
156        elif idType == "pmid":  #FIXME PMID annotation
157            database = self.fromPMID
158        else:
159            return ret
160       
161        for i in ids:
162            if(database.has_key(int(i))):
163                ret[i] = database[int(i)]
164        return ret
165
166
167    def findCIDSubset(self, examples, meshTerms):
168        """
169        Function findCIDSubset will return new list of examples (ids) containing examples which are annotated with at least one MeSH term from \
170        list 'meshTerms'.
171        """
172        newdata = []
173        self.solo_att = self.__findMeshAttribute(examples)
174        ids = list()
175        meshTerms = list(set(meshTerms))
176        # we couldn't find any mesh attribute
177        if self.solo_att == "Unknown":
178            return newdata
179        for e in examples:
180            try:
181                ends = eval(e[self.solo_att].value)
182            except SyntaxError:
183                continue
184            for i in meshTerms:
185                if i in ends:
186                    newdata.append(e['cid'].value)
187                    break
188        return newdata
189   
190    def findFrequentTerms(self, data, minSizeInTerm, treeData=False, callback=None):
191        """
192        Function findFrequentTerms iterates thru examples in data. For every example it finds appropriate annotation (MeSH term). Result of this \
193        function is a dictionary where keys corespond to MeSH terms ids and a value is a number of examples annotated with particular MeSH term.
194        """
195        # we build a dictionary         meshID -> [description, noReference, [cids] ]
196        self.statistics = dict()
197        self.calculated = False
198        self.solo_att = self.__findMeshAttribute(data)
199        # post processing variables
200        ret = dict()
201        ids = []
202        succesors = dict()      # for each term id -> list of succesors
203        succesors["tops"] = []
204        # if we can't identify mesh attribute we return empty data structures
205        if self.solo_att == "Unknown":
206            if treeData:
207                return succesors, ret
208            else:
209                return ret
210        # plain frequency
211        t = 0.0
212        n = len(data)
213        for i in data:
214            t = t + 1
215            if callback and t%10 == 0:
216                callback(int(t*100/n))
217            try:
218                endNodes = list(set(eval(i[self.solo_att].value))) # for every CID we look up for end nodes in mesh. for every end node we have to find its ID 
219            except SyntaxError:
220                #print "Error in parsing ",i[self.solo_att].value
221                continue
222            # we find ID of end nodes
223            endIDs = []
224            for k in endNodes:
225                if(self.toID.has_key(k)):   # this should be always true, but anyway ...
226                    endIDs.extend(self.toID[k])
227                else:
228                    print "Current ontology does not contain MeSH term ", k, "." 
229            # we find id of all parents
230            #allIDs = self.__findParents(endIDs)
231            allIDs = endIDs
232            for k in allIDs:    # for every meshID we update statistics dictionary
233                if(not self.statistics.has_key(k)): # first time meshID
234                    self.statistics[k] = 0
235                self.statistics[k] += 1 # counter increased
236        # post processing
237        for i in self.statistics.iterkeys():
238            if(self.statistics[i] >= minSizeInTerm ): 
239                ret[i] = self.statistics[i]
240                ids.append(i)
241        # we can also return tree data
242        if treeData: #we return also compulsory data for tree printing
243            return self.__treeData(ids), ret
244        else:
245            return ret
246
247    def __treeData(self,ids):
248        succesors = dict()
249        succesors["tops"]= []
250        # we build a list of succesors. for each node we know which are its succesors in mesh ontology
251        for i in ids:
252            succesors[i] = []
253            for j in ids:
254                if(i != j and self.__isPrecedesor(i,j)):
255                    succesors[i].append(j)
256        # for each node from list above we remove its indirect succesors
257        # only  i -1-> j   remain
258        for i in succesors.iterkeys():
259            succs = succesors[i]
260            second_level_succs = []
261            for k in succs: 
262                second_level_succs.extend(succesors[k])
263            for m in second_level_succs:
264                if succesors[i].count(m)>0:
265                    succesors[i].remove(m)
266        # we make a list of top nodes
267        tops = list(ids)
268        for i in ids:
269            for j in succesors[i]:
270                tops.remove(j)
271        # we pack tops table and succesors hash
272        succesors["tops"] = tops
273        return succesors 
274       
275    def findEnrichedTerms(self, reference, cluster, pThreshold=0.05, treeData=False, callback=None, fuzzy=False):
276        """
277        Function findEnrichedTerms computes MeSH term enrichment based on a 'reference' and 'cluster' sets. It returns a dictionary where \
278        keys are enriched (their p-value lower that 'pThreshold') MeSH terms. Key values are lists made of several items (MeSH term id, \
279        MeSH term description, number of examples from the reference set, number of examples from the cluster set, p value, fold enrichment).
280        """
281        self.clu_att = self.__findMeshAttribute(cluster)
282        self.ref_att = self.__findMeshAttribute(reference)
283        if((not self.calculated or self.reference != reference or self.cluster != cluster) and self.ref_att != "Unknown" and \
284        self.clu_att != "Unknown"): # Do have new data? Then we have to recalculate everything.
285            self.reference = reference
286            self.cluster = cluster         
287            self.__calculateAll(callback, fuzzy)
288        # declarations
289        ret = dict()
290        ids = []
291        succesors = dict()      # for each term id -> list of succesors
292        succesors["tops"] = []
293        # if some attributes were unknown
294        if (self.clu_att == "Unknown" or self.ref_att == "Unknown"):
295            if treeData:
296                return  succesors, ret
297            else:
298                return  ret
299        for i in self.statistics.iterkeys():
300            if(self.statistics[i][2] <= pThreshold ) :#or self.statistics[i][4] <= pThreshold ): #
301                ret[i] = self.statistics[i]
302                ids.append(i)
303        if treeData:
304            return self.__treeData(ids),ret
305        else:
306            return ret
307
308    def printMeSH(self, data, selection=["term","r","c", "p"], func=None):
309        """
310        Function printMeSH can be used to print result (dictionary) from the findEnrichedTerms and findFrequentTerms functions.
311        """
312        # first we calculate additional info for printing MeSH ontology
313        info = self.__treeData(data.keys())
314        for i in info["tops"]:
315            self.__pp(0,i,info,data, selection, funct = func)
316
317    def __pp(self, offset, item, relations, data, selection, funct=None):
318        mapping = {"term":0,"desc":1,"r":2,"c":3, "p":4, "fold":5, "func":6} 
319        for i in range(0,offset):
320            print " ",
321        if type(data[item]) == list:
322            pval = "%.4g" % data[item][2]
323            fold = "%.4g" % data[item][3]
324            print_data = [self.toName[item], self.toDesc[self.toName[item]], str(data[item][0]), str(data[item][1]), str(pval), str(fold)]
325            for i in selection:
326                if i != "term":
327                    print i + "=" + print_data[mapping[i]],
328                else:
329                    print print_data[mapping[i]],
330            if funct != None:
331                print " ", funct(print_data[0]),
332            #print self.toName[item], " r=" + str(data[item][1])  +" c="+ str(data[item][2])  ," p=" + str(pval) + " fold=" + str(fold)
333            print ""
334        else:
335            print self.toName[item], " freq=" + str(data[item])
336        for i in relations[item]:
337            self.__pp(offset + 2, i, relations, data, selection, funct = funct) 
338
339    def printHtmlMeSH(self, data, selection=["term","r","c", "p"], func=None):
340        """
341        Function printHtmlMeSH if used to print results (dictinary) from the findFrequentTerms and findFrequentTerms functins in HTML format. \
342        Together with the MeSH ontology it prints data like number of examples, p-values (enrichment).
343        """
344        # first we calculate additional info for printing MeSH ontology
345        info = self.__treeData(data.keys())
346        w = {"term":"'95px'", "r":"'70px'","c":"'70px'","p":"'95px'"}
347        print "<table>\n<tr>"
348        for i in selection:
349            print "<th width=" + w[i] +" align='left'>" + i  +"</th>"
350        if func != None:
351            func("header","")
352        print "</tr>\n"
353        for i in info["tops"]:
354            self.__htmlpp(0,i,info,data, selection, funct = func)
355        print "</table>"
356
357    def __htmlpp(self, offset, item, relations, data, selection, funct = None):
358        mapping = {"term":0,"desc":1,"r":2,"c":3, "p":4, "fold":5, "func":6} 
359        print "<tr>"
360        if type(data[item]) == list:
361            pval = "%.4g" % data[item][2]
362            fold = "%.4g" % data[item][3]
363            print_data = [self.toName[item], self.toDesc[self.toName[item]], str(data[item][0]), str(data[item][1]), str(pval), str(fold)]
364            for i in selection:
365                print "<td>"
366                if i == "term":
367                    for l in range(0,offset):
368                        print "&nbsp;",
369                elif i == "p":
370                    print '%(#)2.3e' % {'#':float(print_data[mapping[i]])} + "</td>",
371                    continue
372                print print_data[mapping[i]] + " &nbsp;</td>",
373
374            if funct != None:
375                print funct(print_data[0],item),
376
377            #print self.toName[item], " r=" + str(data[item][1])  +" c="+ str(data[item][2])  ," p=" + str(pval) + " fold=" + str(fold)
378            print "</tr>"
379        else:
380            print self.toName[item], " freq=" + str(data[item])
381
382        for i in relations[item]:
383            self.__htmlpp(offset + 2, i, relations, data, selection, funct = funct)
384
385    def parsePubMed(self, filename, attributes=["pmid","title","abstract","mesh"], skipExamplesWithout=["mesh"]):
386        """
387        Function parsePubMed can be used to parse (into Orange example table) PubMed search results (in XML).
388        """
389        parser = make_parser()
390        handler = pubMedHandler()
391        parser.setContentHandler(handler)
392        parser.parse(open(filename))
393        atts = []
394        for i in attributes:
395            atts.append(orange.StringVariable(i))
396        domain = orange.Domain(atts,0)
397        data = orange.ExampleTable(domain)
398        print data.domain.attributes
399        mapping = {"pmid":0, "title":1, "abstract":2, "mesh":3, "affilation":4}
400        for i in handler.articles:
401            r=[]
402            skip = False
403            for f in attributes:
404                if skipExamplesWithout.count(f) > 0:
405                    if (f == "mesh" and len(i[mapping[f]]) == 0) or str(i[mapping[f]]) == "":
406                        skip = True
407                r.append(str(i[mapping[f]]))
408            if not skip:
409                data.append(r)
410        return data   
411
412    def __findParents(self,endNodes):
413        """for each end node in endNodes function finds all nodes on the way up to the root"""     
414        res = []
415        for n in endNodes:
416            tmp = n
417            res.append(tmp)
418            for i in range(n.count(".")):
419                tmp = tmp.rstrip("1234567890").rstrip(".")
420                if(tmp not in res):
421                    res.append(tmp)
422        """for i in res:
423            nn = self.toID[self.toName[i]]
424            for t in nn:
425                if not t in res:
426                    res.append(t) """
427        return res
428
429    def __findMeshAttribute(self, data):
430        """
431        Function tries to find attribute which contains list of MeSH terms.
432        """
433        # we get a list of attributes
434        dom = data.domain.attributes
435        for i in dom:             # for each attribute
436            if i.varType == 6:
437                for k in data:         # for each domain
438                    att = str(i.name)
439                    try:                                         # we try to use eval()
440                        r = eval(str(k[att].value))
441                        if type(r) == list:      # attribute type should be list
442                            if self.dataLoaded:      # if ontology is loaded we perform additional test
443                                for i in r:
444                                    if self.toID.has_key(i): return att
445                            else:                  # otherwise we return list attribute
446                                return att
447                    except SyntaxError:
448                        continue
449                    except NameError:
450                        continue   
451        print "Program was unable to determinate MeSH attribute."
452        return "Unknown"
453
454    def __isPrecedesor(self,a,b):
455        """ function returns true if in Mesh ontology exists path from term id a to term id b """
456        if b.count(a) > 0:
457            return True
458        return False
459
460    def __calculateAll(self, callback, fuzzy):
461        """calculates all statistics"""
462        # we build a dictionary         meshID -> [description, noReference,noCluster, enrichment, deprivement, [cids] ]
463        #print "ok"
464        self.statistics = dict()
465        if fuzzy:
466            n = 0
467            for i in self.reference:
468                n = n + float(i['u'])
469            cln = 0
470            for i in self.cluster:
471                cln = cln + float(i['u'])
472
473        else:
474            n = len(self.reference)                                         # reference size
475            cln = len(self.cluster)                                         # cluster size
476        # frequency from reference list
477        r = 0.0
478        for i in self.reference:
479            if callback and r%10 == 0:
480                r += 1.0
481                callback(int(100*r/(n+cln)))
482            try:
483                endNodes = list(set(eval(i[self.ref_att].value))) # for every CID we look up for end nodes in mesh. for every end node we have to find its ID   
484            except SyntaxError:                  # where was a parse error
485                print "Error in parsing ",i[self.ref_att].value
486                if fuzzy:
487                    n = n - float(i["u"])
488                else:
489                    n=n-1
490                continue
491            #we find ID of end nodes
492            endIDs = []
493            for k in endNodes:
494                if(self.toID.has_key(k)):                   # this should be always true, but anyway ...
495                    endIDs.extend(self.toID[k])
496                else:
497                    print "Current ontology does not contain MeSH term ", k, "." 
498            # endIDs may be empty > in this case we can skip this example
499            if len(endIDs) == 0:
500                if fuzzy:
501                    n = n - float(i["u"])
502                else:
503                    n = n-1
504                continue
505            # we find id of all parents
506            #allIDs = self.__findParents(endIDs) PATCH
507            allIDs = endIDs
508            for k in list(set(allIDs)):             # for every meshID we update statistics dictionary
509                if(not self.statistics.has_key(k)):         # first time meshID
510                    self.statistics[k] = [ 0, 0, 0.0, 0.0 ]
511                if fuzzy:
512                    self.statistics[k][0] += float(i["u"])
513                else:
514                    self.statistics[k][0] += 1 # increased noReference
515        # frequency from cluster list
516        r=0.0
517        for i in self.cluster:
518            try:
519                if callback and r%10 == 0:
520                    r += 1.0
521                    callback(int(100*r/(n+cln)))
522                endNodes = list(set(eval(i[self.clu_att].value))) # for every CID we look up for end nodes in mesh. for every end node we have to find its ID   
523            except SyntaxError:
524                #print "Error in parsing ",i[self.clu_att].value
525                if fuzzy:
526                    cln = cln - float(i["u"])
527                else:
528                    cln = cln - 1
529                continue
530            # we find ID of end nodes
531            endIDs = []
532            for k in endNodes:
533                if(self.toID.has_key(k)):               
534                    endIDs.extend(self.toID[k]) # for every endNode we add all corensponding meshIDs
535            # endIDs may be empty > in this case we can skip this example
536            if len(endIDs) == 0:
537                if fuzzy:
538                    cln = cln - float(i["u"])
539                else:
540                    cln = cln-1
541                continue
542            # we find id of all parents
543            #allIDs = self.__findParents(endIDs)
544            allIDs = endIDs
545            for k in list(set(allIDs)): # for every meshID we update statistics dictionary
546                if self.statistics.has_key(k):
547                    if fuzzy:
548                        self.statistics[k][1] += float(i["u"])
549                    else:
550                        self.statistics[k][1] += 1 # increased noCluster
551        self.ratio = float(cln)/float(n)
552        # enrichment
553        for i in self.statistics.iterkeys():
554            self.statistics[i][2] = self.__calcEnrichment(int(n),int(cln),int(self.statistics[i][0]),int(self.statistics[i][1]))    # p enrichment
555            self.statistics[i][3] = float(self.statistics[i][1]) / float(self.statistics[i][0] ) / self.ratio   # fold enrichment
556        self.calculated = True     
557
558    def __calcEnrichment(self, n, c, t, tc):
559        """
560        Function calculates hypergeometric distribution based on the input parameters.
561        n - total number of items ie. size(cluster + reference)
562        c - cluster size ie. size(cluster)
563        t - number of all items in observed MeST term group
564        tc - number of cluster chemicals in observed term group
565        """
566
567        # FIXME: Popravi cudno racunanje enrichmenta v mejnih primerih.
568       
569        result=0
570        for i in range(0,tc):
571            result = result + exp(self.log_comb(t,i) + self.log_comb(n-t,c-i))
572        result = result*1.0 / exp(self.log_comb(n,c))
573        return (1.0-result)
574
575    def log_comb(self,n, m): # it returns log(nCr(n,m))
576        return self.lookup[n] - self.lookup[n-m] - self.lookup[m]
577
578    def __loadOntologyFromDisk(self):
579        """
580        Function loads MeSH ontology and chemical annotation into internal data structures.
581        """
582        self.toID = dict()  #   name -> [IDs] Be careful !!! One name can match many IDs!
583        self.toName = dict() #   ID -> name
584        self.toDesc = dict() #  name -> description
585        self.fromCID = dict() #   cid -> term id
586        self.fromPMID = dict()  #   pmid -> term id
587       
588        __dataPath = orngServerFiles.localpath('MeSH')
589        try:       
590            # reading graph structure from file
591            d = file(os.path.join(__dataPath,'mesh-ontology.dat'))
592        except IOError:
593            print os.path.join(__dataPath,'mesh-ontology.dat') + " does not exist!"
594            return False
595        try:   
596            # reading cid annotation from file
597            f = file(os.path.join(__dataPath,'cid-annotation.dat'))
598        except IOError:
599            print os.path.join(__dataPath,'cid-annotation.dat') + " does not exist!"
600            #return False
601        # loading ontology graph
602        t=0
603        for i in d:
604            t += 1
605            parts = i.split("\t") # delimiters are tabs
606            if(len(parts) != 3):
607                print "error reading ontology ", parts[0]
608            parts[2] = parts[2].rstrip("\n\r")
609            ids = parts[1].split(";")
610            self.toID[parts[0]] = ids   # append additional ID
611            self.toDesc[parts[0]] = parts[2]
612            for r in ids:
613                self.toName[r] = parts[0]
614            # loading cid -> mesh
615            for i in f:
616                parts = i.split(";")        # delimiters are tabs
617                if(len(parts) != 2):
618                    print "error reading ontology ", parts[0]
619                parts[1] = parts[1].rstrip("\n\r")
620                cid = int(parts[0])
621                if self.fromCID.has_key(cid):
622                    self.fromCID[cid].append(parts[1])
623                else:
624                    self.fromCID[cid] = [parts[1]]
625            # loading pmid -> mesh, TODO
626        #print "Current MeSH ontology contains ", t, " mesh terms."
627        return True
628
629class pubMedHandler(ContentHandler):
630    def __init__(self):
631        self.state = 0  #  0 start state, 1 pmid, 2 title, 3 abstract, 4 mesh
632        self.articles = []
633        self.pmid = "0"
634        self.title = ""
635        self.mesh = list()
636        self.abstract = ""
637        self.affiliation = ""
638
639    def startElement(self, name, attributes):
640        # print "parsam ", name
641        if name == "PubmedArticle":
642            self.pmid = ""
643            self.abstract = ""
644            self.title = ""
645            self.mesh = []
646        if name == "PMID":
647            self.state = 1
648        if name == "ArticleTitle":
649            self.state = 2
650        if name == "AbstractText":
651            self.state = 3
652        if name == "DescriptorName":
653            self.state = 4
654            self.mesh.append("")
655        if name == "Affiliation":
656            self.state = 5
657
658    def characters(self, data):
659        if self.state == 1:
660            self.pmid += data
661        if self.state == 2:
662            self.title += data.encode("utf-8")
663        if self.state == 3:
664            self.abstract += data.encode("utf-8")
665        if self.state == 4:
666            self.mesh[-1] += data.encode("utf-8")
667        if self.state == 5:
668            self.affiliation += data.encode("utf-8")
669
670    def endElement(self, name):
671        #print "   koncujem ", name
672        self.state = 0
673        if name == "PubmedArticle":
674            self.articles.append([self.pmid,self.title,self.abstract,self.mesh, self.affiliation])
675
676class MappedMeSHParser(SGMLParser):
677    def reset(self):
678        self.pieces = []
679        self.terms = []
680        self.foundMeSH = False
681        self.nextIsTerm = False
682        self.endTags = ['Display','Write to the Help Desk', 'Previous Indexing:', 'Entry Terms:','Pharmacologic Action:' ]
683        SGMLParser.reset(self)
684   
685    def unknown_starttag(self, tag, attrs):
686        strattrs = "".join([' %s="%s"' % (key, value) for key, value in attrs])
687        if self.foundMeSH and tag=='a':
688            self.nextIsTerm = True
689   
690    def handle_data(self, text):
691        text = text.strip()
692        if text == '':
693            return
694        if text == 'Heading Mapped to:':
695            self.foundMeSH = True
696        if self.endTags.count(text)>0:
697            self.foundMeSH = False
698        elif self.nextIsTerm:
699            self.terms.append(text)
700            self.nextIsTerm = False
701
702class PubChemMeSHParser(SGMLParser):
703    def reset(self):
704        self.next = 0
705        self.nextLink = ''
706        self.directTerms = []
707        self.indirectTerms = []
708        self.foundMeSH = False
709        self.foundIndirectMeSH = 0
710        SGMLParser.reset(self)
711        # strategy as follows
712        # Beetween strings "Drug and Chemical Info" and ("Pharmalogical Action" or "PubMed via MeSH" or "PubMed MeSH Keyword Summary") find hyperlinks. Based on title attribute we can distingue direct MeSH terms beetween mapped terms.
713
714    def unknown_starttag(self, tag, attrs): 
715        #if self.foundMeSH:
716        #print 'tag ', tag, attrs
717        if tag=='a' and len(attrs) > 2 and attrs[0][0] == 'name' and attrs[0][1] == 'gomesh': #print attrs
718            self.nextLink = attrs[1][1]
719
720    def handle_data(self, text):
721        text = text.strip()
722        if text == '':
723            return
724        #print 'data ', text
725        if self.foundMeSH:
726            #print '*'
727            self.directTerms.append(text)
728       
729        if self.foundIndirectMeSH == 2:
730            #print "indirect term", text, self.nextLink
731            self.indirectTerms.append((text,self.nextLink))
732            self.foundIndirectMeSH = 0
733       
734        if self.foundIndirectMeSH == 1:
735            self.foundIndirectMeSH = 2
736               
737        if text == 'Drug and Chemical Information:':
738            self.foundIndirectMeSH = 1
739
740
741        if text == 'Pharmacological Action' or text == 'Medication Information':
742            self.foundIndirectMeSH = 0
743        if text == "Chemical Classification":
744            self.foundMeSH = True
745        elif (text == "Safety and Toxicology" or text=="Literature" or text == "Classification" or text == "PubMed via MeSH" or text == "PubMed MeSH Keyword Summary"):
746            self.foundMeSH = False
747            if len(self.directTerms) > 0:
748                self.directTerms.pop()
749
750class SmilesParser(SGMLParser):
751    def reset(self):
752        self.nextSmile = False
753        self.smiles = ""
754        SGMLParser.reset(self)
755
756    def unknown_starttag(self, tag, attrs): 
757        #print tag, " ", attrs
758        if tag=="item":
759            for (i,j) in attrs:
760                if i=="name" and j=="CanonicalSmile":
761                    self.nextSmile = True
762
763    def handle_data(self, text):
764        if self.nextSmile:
765            self.smiles = text
766            self.nextSmile = False
767
768class CIDsParser(SGMLParser):
769    def reset(self):
770        self.cid = False
771        self.cids = []
772        SGMLParser.reset(self)
773
774    def unknown_starttag(self, tag, attrs): 
775        #print tag, " ", attrs
776        if tag=="id":
777            self.cid = True
778           
779    def unknown_endtag(self, tag): 
780        #print tag,
781        if tag=="id":
782            self.cid = False
783
784    def handle_data(self, text):
785        #print text
786        if self.cid and not text == '\n\t\t':
787            #print int(text.strip())
788            self.cids.append(int(text))
789
790class FormulaParser(SGMLParser):
791    def reset(self):
792        self.nextFormula = False
793        self.formula = ""
794        SGMLParser.reset(self)
795
796    def unknown_starttag(self, tag, attrs): 
797        #print tag, " ", attrs
798        if tag=="item":
799            for (i,j) in attrs:
800                if i=="name" and j=="MolecularFormula":
801                    self.nextFormula = True
802
803    def handle_data(self, text):
804        if self.nextFormula:
805            self.formula = text
806            self.nextFormula = False
807
808class CIDSParser(SGMLParser):
809    def reset(self):
810        self.nextCID = False
811        self.cids = []
812        SGMLParser.reset(self)
813
814    def unknown_starttag(self, tag, attrs): 
815        #print tag, " ", attrs
816        if tag=="item":
817            for (i,j) in attrs:
818                if i=="name" and j=="int":
819                    self.nextCID = True
820
821    def handle_data(self, text):
822        if self.nextCID:
823            self.cids.append(text)
824            self.nextCID = False
825
826class pubChemAPI(object):
827    """
828    Class pubChemAPI is used to simplift the access to the pubChem database.
829    """
830    def __init__(self):
831        self.data = ""
832        self.identifier = ""
833        self.database = ""
834
835    def getSMILE(self, id, typ):
836        """
837        Function getSMILE takes chemicals id and its type and returns a SMILES representation.
838        """
839        dbs = {"sid":"pcsubstance", "cid":"pccompound"}
840        if not dbs.has_key(typ):
841            return "Unknown identifier"
842        # maybe we already have the data ...
843        if id == self.identifier and dbs[typ] == self.database:
844            raw = self.data
845        else:
846            url2fetch = "http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=" + dbs[typ] + "&id=" + str(id)
847            #print url2fetch
848            d = urlopen(url2fetch)
849            raw = d.read()
850            self.data = raw
851            self.identifier = id
852            self.database = dbs[typ]
853        pr = SmilesParser()
854        pr.feed(raw)
855        data = pr.smiles
856        pr.close()
857        return data
858
859    def getMeSHterms(self, cid):   
860        """
861        Functions getMeSHterms tries to find MeSH annotation for given the CID in the pubChem database.
862        """
863        usock = urlopen("http://pubchem.ncbi.nlm.nih.gov/summary/summary.cgi?cid="+ str(cid)  +"&viewopt=PubChem&hcount=100&tree=t#MeSH")
864        parser = PubChemMeSHParser()
865        text = usock.read()
866        #print text
867        parser.feed(text)
868        usock.close()
869        parser.close()
870        allTerms = parser.directTerms
871        #print allTerms
872        for (k,i) in parser.indirectTerms:
873            #print k, i
874            #continue
875            usock = urlopen(i)
876            parser = MappedMeSHParser()
877            parser.feed(usock.read())
878            allTerms.extend(parser.terms)   
879            usock.close()
880            parser.close()
881        # from allTerms make a set
882        allTerms = list(set(allTerms))
883        return allTerms
884
885    def getMolecularFormula(self, id, typ):
886        """
887        Functions getMolecularFormula tries to find molecular formula for the given the chemical id and type in the pubChem database.
888        """
889        dbs = {"sid":"pcsubstance", "cid":"pccompound"}
890        if not dbs.has_key(typ):
891            return "Unknown identifier"
892        # maybe we already have the data ...
893        if id == self.identifier and dbs[typ] == self.database:
894            raw = self.data
895        else:
896            url2fetch = "http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=" + dbs[typ] + "&id=" + str(id)
897            #print url2fetch
898            d = urlopen(url2fetch)
899            raw = d.read()
900            self.data = raw
901            self.identifier = id
902            self.database = dbs[typ]
903        pr = FormulaParser()
904        pr.feed(raw)
905        data = pr.formula
906        pr.close()
907        return data
908
909    def getCIDfromSID(self, sid):
910        """
911        Functions getCIDfromSID converts compound id to substance id.
912        """
913        dbs = {"sid":"pcsubstance", "cid":"pccompound"}
914        if not dbs.has_key("sid"):
915            return "Unknown identifier"
916        # maybe we already have the data ...
917        if sid == self.identifier and dbs["sid"] == self.database:
918            raw = self.data
919        else:
920            url2fetch = "http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=" + dbs["sid"] + "&id=" + str(sid)
921            #print url2fetch
922            d = urlopen(url2fetch)
923            raw = d.read()
924            self.data = raw
925            self.identifier = sid
926            self.database = dbs["sid"]
927        pr = CIDSParser()
928        pr.feed(raw)
929        data = pr.cids
930        pr.close()
931        return data
932
933    def getPharmActionList(self, cid):
934        """
935        Functions getPharmActionList finds pharmacological MeSH annotation for the given CID.
936        <Item Name="PharmActionList" Type="List">
937        """ 
938        dbs = {"sid":"pcsubstance", "cid":"pccompound"}
939        # maybe we already have the data ...
940        if cid == self.identifier and dbs["cid"] == self.database:
941            raw = self.data
942        else:
943            url2fetch = "http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=" + dbs["cid"] + "&id=" + str(cid)
944            d = urlopen(url2fetch)
945            raw = d.read()
946        #print url2fetch
947        self.data = raw
948        self.identifier = cid
949        self.database = dbs["cid"]
950        ndata = raw.replace('\t','').replace('\n','')
951        xmldoc = minidom.parseString(ndata)
952        items = xmldoc.getElementsByTagName('Item')
953        data = []
954        for i in items:
955            if i.attributes['Name'].value == 'PharmActionList':
956                terms = i.childNodes
957                for t in terms:
958                    data.append(str(t.childNodes[0].data))
959        return data
960
961    def getCIDs(self, name, weight):
962        """
963        Functions getCIDs tries to find correct chemical id based on given chemical name and weight.
964        """
965        offset = 0.005
966        url2fetch = "http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=pccompound&term="
967        url2fetch += name.replace(" ","%20") + "%20" + str(weight*(1-offset)) + ":" + str(weight*(1+offset)) + "[mw]&retmax=10000"
968        #print url2fetch
969        d = urlopen(url2fetch)
970        raw = d.read()
971        e = CIDsParser()
972        e.feed(raw)
973        return e.cids
974
975    def getCIDs(self,name):
976        """
977        Functions getCIDs tries to find correct chemical id based on given chemical name.
978        """
979        #print "k"
980        url2fetch = "http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=pccompound&term="
981        url2fetch += name.replace(" ","%20") + "&retmax=10000"
982        #print url2fetch
983        d = urlopen(url2fetch)
984        raw = d.read()
985        e = CIDsParser()
986        e.feed(raw)
987        return e.cids
988
989    def __strip_ml_tags__(self,in_text):
990        s_list = list(in_text)
991        i,j = 0,0
992        while i < len(s_list):
993            if s_list[i] == '<':
994                while s_list[i] != '>':
995                    s_list.pop(i)
996                s_list.pop(i)
997            else:
998                i=i+1
999        join_char=''
1000        return join_char.join(s_list)
Note: See TracBrowser for help on using the repository browser.