source: orange-bioinformatics/Orange/bioinformatics/obiMeSH.py @ 1625:cefeb35cbfc9

Revision 1625:cefeb35cbfc9, 30.8 KB checked in by mitar, 2 years ago (diff)

Moving files around.

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