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

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

Moving files around.

Line 
1import urllib2
2import os, sys, shutil
3import cPickle
4import tarfile
5import StringIO
6import obiGenomicsUpdate
7import obiData
8import orngEnviron
9import orngServerFiles
10
11from urllib2 import urlopen
12from collections import defaultdict
13
14# list of common organisms from http://www.ncbi.nlm.nih.gov/Taxonomy
15def common_taxids():
16    """Return taxonomy IDs for common organisms."""
17    return ["3702",  # Arabidopsis thaliana
18            "9913",  # Bos taurus
19            "6239",  # Caenorhabditis elegans
20            "3055",  # Chlamydomonas reinhardtii
21            "7955",  # Danio rerio (zebrafish)
22            "352472", # Dictyostelium discoideum
23            "7227",  # Drosophila melanogaster
24            "562",   # Escherichia coli
25            "11103", # Hepatitis C virus
26            "9606",  # Homo sapiens
27            "10090", # Mus musculus
28            "2104",  # Mycoplasma pneumoniae
29            "4530",  # Oryza sativa
30            "5833",  # Plasmodium falciparum
31            "4754",  # Pneumocystis carinii
32            "10116", # Rattus norvegicus
33            "4932",  # Saccharomyces cerevisiae
34            "4896",  # Schizosaccharomyces pombe
35            "31033", # Takifugu rubripes
36            "8355",  # Xenopus laevis
37            "4577",  # Zea mays
38             ] 
39
40def essential_taxids():
41    """Return taxonomy IDs for organisms that are included in (default) Orange Bioinformatics installation."""
42    return ["352472", # Dictyostelium discoideum
43            "7227",  # Drosophila melanogaster
44            "9606",  # Homo sapiens
45            "10090", # Mus musculus
46            "4932",  # Saccharomyces cerevisiae
47            ] 
48
49default_database_path = os.path.join(orngEnviron.bufferDir, "bigfiles", "Taxonomy")
50
51class MultipleSpeciesException(Exception):
52    pass
53
54class UnknownSpeciesIdentifier(Exception):
55    pass
56
57def pickled_cache(filename=None, dependencies=[], version=1, maxSize=30):
58    """ Return a cache function decorator.
59    """
60    def datetime_info(domain, filename):
61        try:
62            return orngServerFiles.info(domain, filename)["datetime"]
63        except IOError:
64            return orngServerFiles.ServerFiles().info(domain, filename)["datetime"]
65           
66    def cached(func):
67        default_filename = os.path.join(orngEnviron.bufferDir, func.__module__ + "_" + func.__name__ + "_cache.pickle")
68        def f(*args, **kwargs):
69            currentVersion = tuple([datetime_info(domain, file) for domain, file in dependencies]) + (version,)
70            try:
71                cachedVersion, cache = cPickle.load(open(filename or default_filename, "rb"))
72                if cachedVersion != currentVersion:
73                    cache = {}
74            except IOError, er:
75                cacheVersion, cache = "no version", {}
76            allArgs = args + tuple([(key, tuple(value) if type(value) in [set, list] else value)\
77                                     for key, value in kwargs.items()])
78            if allArgs in cache:
79                return cache[allArgs]
80            else:
81                res = func(*args, **kwargs)
82                if len(cache) > maxSize:
83                    del cache[iter(cache).next()]
84                cache[allArgs] = res
85                cPickle.dump((currentVersion, cache), open(filename or default_filename, "wb"), protocol=cPickle.HIGHEST_PROTOCOL)
86                return res
87        return f
88
89    return cached
90
91def cached(func):
92    """Cached one arg method
93    """
94    def f(self, arg):
95        if arg not in self._cache:
96            self._cache[arg] = func(self, arg)
97        return self._cache[arg]
98    f._cache = {}
99    f.__name__ = "Cached " + func.__name__
100    return f
101   
102class TextDB(object):
103    entry_start_string = chr(255)
104    entry_end_string = chr(254)+"\n"
105    entry_separator_string = chr(253)
106   
107    @property
108    def _text_lower(self):
109        if id(self._text) == self._lower_text_id:
110            return self._lower_text
111        else:
112            self._lower_text_id = id(self._text)
113            self._lower_text = self._text.lower()
114            return self._lower_text
115       
116    def __init__(self, file=None, **kwargs):
117        self._text = ""
118        self._lower_text_id = id(self._text) - 1
119        self._cache = {}
120        self.__dict__.update(kwargs)
121       
122        if file != None:
123            self._text = open(file, "rb").read()
124
125    def _find_all(self, string, start=0, text=None, unique=True):
126        text = text if text != None else self._text_lower
127        while True:
128            index = text.find(string, start)
129            if index != -1:
130                yield index
131                if unique:
132                    start = text.find(self.entry_start_string, index + 1)
133                else:
134                    start = index + 1
135            else:
136                raise StopIteration
137
138    def _get_entry_at(self, index, text=None):
139        text = text if text != None else self._text
140        start = text.rfind(self.entry_start_string, 0, index + 1)
141        end = text.find(self.entry_end_string, index)
142        return self._text[start+1:end]
143
144    @cached
145    def get_entry(self, id):
146        try:
147            index = self._find_all(self.entry_start_string + id + self.entry_separator_string).next()
148        except StopIteration:
149            raise KeyError, id
150        return self._get_entry_at(index)
151               
152    def search(self, string):
153        string = string.lower()
154        res = []
155        for idx in self._find_all(string):
156            entry = self._get_entry_at(idx)
157            id , rest = entry.split(self.entry_separator_string, 1)
158            self._cache[id] = entry
159            res.append(id)
160        return res
161
162    def insert(self, entry):
163        self._text += self.entry_start_string + self.entry_separator_string.join(entry) + self.entry_end_string
164
165    def __iter__(self):
166        for idx in self._find_all(self.entry_start_string):
167            entry = self._get_entry_at(idx)
168            if entry:
169                yield entry.split(self.entry_separator_string ,1)[0]
170
171    def __getitem__(self, id):
172        entry = self.get_entry(id)
173        return entry.split(self.entry_separator_string)[1:]
174
175    def __setitem__(self, id, entry):
176        self.insert([id] + list(entry))
177
178    def create(self, filename):
179        f = open(filename, "wb")
180        f.write(self._text)
181        def write(entry):
182            f.write(self.entry_start_string + self.entry_separator_string.join(entry) + self.entry_end_string)
183        return write
184   
185class Taxonomy(object):
186    __shared_state = {"_text":None, "_info":None}
187    def __init__(self):
188        self.__dict__ = self.__shared_state
189        if not self._text:
190            self.Load()
191           
192    def Load(self):
193        try:
194            self._text = TextDB(os.path.join(default_database_path, "ncbi_taxonomy.tar.gz", "ncbi_taxonomy.db"))
195            self._info = TextDB(os.path.join(default_database_path, "ncbi_taxonomy.tar.gz", "ncbi_taxonomy_inf.db"))
196            return
197        except Exception, ex:
198            pass
199        try:
200            import orngServerFiles as sf
201            sf.download("Taxonomy", "ncbi_taxonomy.tar.gz")
202            self._text = TextDB(os.path.join(default_database_path, "ncbi_taxonomy.tar.gz", "ncbi_taxonomy.db"))
203            self._info = TextDB(os.path.join(default_database_path, "ncbi_taxonomy.tar.gz", "ncbi_taxonomy_inf.db"))
204            return
205        except Exception, ex:
206            raise
207
208    def get_entry(self, id):
209        try:
210            entry = self._text[id]
211        except KeyError:
212            raise UnknownSpeciesIdentifier
213        return entry
214               
215    def search(self, string, onlySpecies=True):
216        res = self._text.search(string)
217        if onlySpecies:
218            res = [r for r in res if "species" in self._text[r][1]]
219        return res
220
221    def __iter__(self):
222        return iter(self._text)
223
224    def __getitem__(self, id):
225        entry = self.get_entry(id)
226        return entry[2] ## item with index 2 is allways scientific name
227
228    def other_names(self, id):
229        entry = self.get_entry(id)
230        info = self._info[id]
231        names = entry[2:] ## index 2 and larger are names
232        return list(zip(names, info))[1:] ## exclude scientific name
233
234    def rank(self, id):
235        entry = self.get_entry(id)
236        return entry[1]
237
238    def parent(self, id):
239        entry = self.get_entry(id)
240        return entry[0]
241
242    def subnodes(self, id, levels=1):
243        res = self._text.search(self._text.entry_separator_string + id + self._text.entry_separator_string)
244        res = [r for r in res if self.get_entry(r)[0] == id]
245        if levels > 1:
246            for r in list(res):
247                res.extend(self.subnodes(r, levels-1))
248        return res
249
250    def taxids(self):
251        return list(self)
252   
253    @staticmethod
254    def ParseTaxdumpFile(file=None, outputdir=None, callback=None):
255        from cStringIO import StringIO
256        if file == None:
257            file = tarfile.open(None, "r:gz", StringIO(urlopen("ftp://ftp.ncbi.nih.gov/pub/taxonomy/taxdump.tar.gz").read()))
258        if type(file) == str:
259            file = tarfile.open(file)
260        names = file.extractfile("names.dmp").readlines()
261        nodes = file.extractfile("nodes.dmp").readlines()
262        namesDict = defaultdict(list)
263        for line in names:
264            if not line.strip():
265                continue
266            line = line.rstrip("\t\n|").split("\t|\t")
267            id, name, unique_name, name_class = line
268            if unique_name:
269                namesDict[id].append((unique_name , name_class))
270            else:
271                namesDict[id].append((name , name_class))
272
273        nodesDict = {}
274        for line in nodes:
275            if not line.strip():
276                continue
277            line = line.split("\t|\t")[:3]
278            id, parent, rank = line
279            nodesDict[id] = (parent, rank)
280       
281        name_class_codes = defaultdict(iter(range(255)).next)
282        name_class_codes["scientific name"] ## Force scientific name to be first
283        if outputdir == None:
284            outputdir = default_database_path
285        text = TextDB().create(os.path.join(outputdir, "ncbi_taxonomy.db"))
286        info = TextDB().create(os.path.join(outputdir, "ncbi_taxonomy_inf.db"))
287        milestones = set(range(0, len(namesDict), max(len(namesDict)/100, 1)))
288        for i, (id, names) in enumerate(namesDict.items()):
289            parent, rank = nodesDict[id]
290            ## id, parent and rank go first
291            entry = [id, parent, rank]
292            ## all names and name class codes pairs follow ordered so scientific name is first
293            names = sorted(names, key=lambda (name, class_): name_class_codes[class_])
294            entry.extend([name for name ,class_ in names])
295            info_entry = [id] + [class_ for name, class_ in names]
296            text(entry)
297            info(info_entry)
298            if callback and i in milestones:
299                callback(i)
300   
301def name(taxid):
302    """ Return the scientific name for organism with taxid.
303    """
304    return Taxonomy()[taxid]
305
306def other_names(taxid):
307    """ Return a list of (name, name_type) tuples but exclude the scientific name.
308    """
309    return  Taxonomy().other_names(taxid)
310
311@pickled_cache(None, [("Taxonomy", "ncbi_taxonomy.tar.gz")], version=1)
312def search(string, onlySpecies=True, exact=False):
313    """ Search the NCBI taxonomy database for an organism
314    Arguments::
315            - *string*      Search string
316            - *onlySpecies* Return only taxids of species (and subspecies)
317            - *exact*       Return only taxids of organism that exactly match the string
318    """
319    ids = Taxonomy().search(string, onlySpecies)
320    if exact:
321        ids = [id for id in ids if string in [name(id)] + [t[0] for t in other_names(id)]]
322    return ids
323
324def lineage(taxid):
325    """ Return a list of taxids ordered from the topmost node (root) to taxid.
326    """
327    tax = Taxonomy()
328    result = [taxid]
329    while True:
330        parent = tax.parent(result[-1])
331        result.append(parent)
332        if tax[parent] == "root" or parent=="1":
333            break
334    result.reverse()
335    return result
336   
337def to_taxid(code, mapTo=None):
338    """ See if the code is a valid code in any database and return a set of its taxids.
339    """
340    import obiKEGG, obiGO
341    results = set()
342    for test in [obiKEGG.to_taxid, obiGO.to_taxid]:
343        try:
344            r = test(code)
345            if type(r) == set:
346                results.update(r)
347            else:
348                results.add(r)
349        except Exception, ex:
350            pass
351
352    if mapTo:
353        mapped = [[parent_id for parent_id in mapTo if parent_id in lineage(r)].pop() for r in results if any(parent_id in lineage(r) for parent_id in mapTo)]
354        if not mapped:
355            t = Taxonomy()
356            subnodes = dict([(r, t.subnodes(r, 4)) for r in results])
357            mapped = [id for id in mapTo if any(id in subnodes[r] for r in results)]
358        results = mapped
359
360    return results
361
362def taxids():
363    return Taxonomy().taxids()
364
365import obiGenomicsUpdate
366
367class Update(obiGenomicsUpdate.Update):
368    def GetDownloadable(self):
369        return [Update.UpdateTaxonomy]
370   
371    def IsUpdatable(self, func, args):
372        from datetime import datetime
373        import obiData
374        if func == Update.UpdateTaxonomy:
375##            stream = urllib2.urlopen("ftp://ftp.ncbi.nih.gov/pub/taxonomy/taxdump.tar.gz")
376##            date = datetime.strptime(stream.headers.get("Last-Modified"), "%a, %d %b %Y %H:%M:%S %Z")
377            ftp = obiData.FtpWorker("ftp.ncbi.nih.gov")
378            size, date = ftp.statFtp("pub/taxonomy/taxdump.tar.gz")
379            return date > self.GetLastUpdateTime(func, args)
380
381    def UpdateTaxonomy(self):
382        Taxonomy.ParseTaxdumpFile(outputdir=self.local_database_path)
383        import tarfile
384        tFile = tarfile.open(os.path.join(self.local_database_path, "ncbi_taxonomy.tar.gz"), "w:gz")
385        tFile.add(os.path.join(self.local_database_path, "ncbi_taxonomy.db"), "ncbi_taxonomy.db")
386        tFile.add(os.path.join(self.local_database_path, "ncbi_taxonomy_inf.db"), "ncbi_taxonomy_inf.db")
387        tFile.close()
388
389if __name__ == "__main__":
390    ids = search("Homo sapiens")
391    print ids
392    print other_names(ids[0])
393   
Note: See TracBrowser for help on using the repository browser.