source: orange-bioinformatics/Orange/bioinformatics/obiMeSH.py @ 1632:9cf919d0f343

Revision 1632:9cf919d0f343, 30.8 KB checked in by mitar, 2 years ago (diff)

Fixing imports.

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