Changeset 1504:ce65c16ede61 in orange-bioinformatics


Ignore:
Timestamp:
10/27/11 18:04:20 (2 years ago)
Author:
ales_erjavec <ales.erjavec@…>
Branch:
default
Convert:
5aacf9cbea4421aa1fdc56b11f54a56bc7c171d5
Message:

Added PPIDatabase like interface for GeneMania.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • obiGeneMania.py

    r1288 r1504  
    389389    return Connection().retrieve(org, genes, m, r) 
    390390 
     391 
     392""" 
     393====================== 
     394PPI Database interface 
     395====================== 
     396 
     397""" 
     398 
     399import obiPPI 
     400import orngServerFiles 
     401import sqlite3 
     402import csv 
     403import os 
     404import posixpath 
     405 
     406from contextlib import contextmanager 
     407import StringIO 
     408 
     409@contextmanager 
     410def finishing(obj): 
     411    """ Calls obj.finish() on context exit. 
     412    """ 
     413    yield obj 
     414    obj.finish() 
     415 
     416def guess_size(fileobj): 
     417    try: 
     418        if isinstance(fileobj, file): 
     419            return os.fstat(fileobj.fileno()).st_size 
     420        elif isinstance(fileobj, StringIO.StringIO): 
     421            pos = fileobj.tell() 
     422            fileobj.seek(0, 2) 
     423            length = fileobj.tell() - pos 
     424            fileobj.seek(pos, 0) 
     425            return length 
     426        elif isinstance(fileobj, urllib.addinfourl): 
     427            length = fileobj.headers.get("content-length", None) 
     428            return length 
     429    except Exception, ex: 
     430        pass 
     431 
     432 
     433def copyfileobj(src, dst, buffer=2**10, content_len=None, progress=None): 
     434    count = 0 
     435    if content_len is None: 
     436        content_len = guess_size(src) or sys.maxint 
     437    while True: 
     438        data = src.read(buffer) 
     439        dst.write(data) 
     440        count += len(data) 
     441        if progress: 
     442            progress(100.0 * count / content_len) 
     443        if not data: 
     444            break 
     445             
     446             
     447def wget(url, directory=".", dst_obj=None, progress=None): 
     448    """ 
     449    .. todo:: Move too Orange.misc 
     450    """ 
     451    stream = urllib2.urlopen(url) 
     452    length = stream.headers.get("content-length", None) 
     453    if length is None: 
     454        length = sys.maxint 
     455    else: 
     456        length = int(length) 
     457     
     458    basename = posixpath.basename(url) 
     459         
     460    if dst_obj is None: 
     461        dst_obj = open(os.path.join(directory, basename), "wb") 
     462     
     463    if progress == True: 
     464        from Orange.misc import ConsoleProgressBar 
     465        progress = ConsoleProgressBar("Downloading %r." % basename) 
     466        with finishing(progress): 
     467            copyfileobj(stream, dst_obj, buffer=2**10, content_len=length, 
     468                        progress=progress) 
     469    else: 
     470        copyfileobj(stream, dst_obj, buffer=2**10, content_len=length, 
     471                    progress=progress) 
     472     
     473import obiTaxonomy 
     474from collections import namedtuple 
     475from operator import itemgetter 
     476from Orange.misc import lru_cache 
     477 
     478GENE_MANIA_INTERACTION_FIELDS = \ 
     479    ["gene_a", "gene_b", "weight", "network_name", 
     480     "network_group", "source", "pubmed_id"] 
     481      
     482GeneManiaInteraction = namedtuple("GeneManiaInteraction", 
     483                                  field_names=GENE_MANIA_INTERACTION_FIELDS 
     484                                 ) 
     485 
     486import weakref 
     487class Internalizer(object): 
     488    """ A class that acts as the python builtin function ``intern``, 
     489    for as long as it is alive. 
     490     
     491    .. note:: This is for memory optimization only, it does not affect  
     492        dict lookup speed. 
     493     
     494    """ 
     495    def __init__(self): 
     496        self._intern_dict = {} 
     497         
     498    def __call__(self, obj): 
     499        return self._intern_dict.setdefault(obj, obj) 
     500     
     501class GeneManiaDatabase(obiPPI.PPIDatabase): 
     502    DOMAIN = "PPI" 
     503    SERVER_FILE = "gene-mania-{taxid}.sqlite" 
     504     
     505    TAXID2NAME = "" 
     506     
     507    # DB schema 
     508    SCHEMA = """ 
     509    table: `genes` 
     510        - `internal_id`: int (pk) 
     511        - `gene_name`: text (preferred name) 
     512         
     513    table: `synonyms`: 
     514        - `internal_id: int (foreign key `genes.internal_id`) 
     515        - `synonym`: text 
     516        - `source_id`: int 
     517         
     518    table: `source`: 
     519        - `source_id`: int 
     520        - `source_name`: text 
     521         
     522    table: `links`: 
     523        - `gene_a`: int (fk `genes.internal_key`) 
     524        - `gene_b`: int (fk `genes.internal_key`) 
     525        - `network_id`: (fk `networks.network_id`) 
     526        - `weight`: real 
     527         
     528    table: `networks`: 
     529        - `network_id`: int 
     530        - `network_name`: text 
     531        - `network_group`: text 
     532        - `source`: text 
     533        - `pubmed_id`: text 
     534         
     535    view: `links_annotated`: 
     536        - `gene_name_a` 
     537        - `gene_name_b` 
     538        - `network_name` 
     539        - `network_group` 
     540        - `weight` 
     541         
     542    """ 
     543     
     544    def __init__(self, taxid): 
     545        self.taxid = taxid 
     546         
     547    @classmethod 
     548    def common_taxids(self): 
     549        return ["3702", "6239", "7227", "9606", "10090", "10116", "4932"] 
     550     
     551    def organisms(self): 
     552        """ Return all organism taxids contained in this database. 
     553         
     554        .. note:: a single taxid is returned (the one used at 
     555            instance initialization)    
     556         
     557        """ 
     558        return [self.taxid] 
     559     
     560    def ids(self, taxid=None): 
     561        """ Return all primary ids for `taxid`. 
     562        """ 
     563        if taxid is None: 
     564            taxids = self.organisms() 
     565            return reduce(list.__add__, map(self.ids, taxids), []) 
     566         
     567        con = self._db(taxid) 
     568        cur = con.execute("""\ 
     569            SELECT gene_name FROM genes 
     570            """) 
     571        return map(itemgetter(0), cur) 
     572         
     573    def synonyms(self, id): 
     574        """ Return a list of synonyms for primary `id`. 
     575        """ 
     576        con = self._db(self.taxid) 
     577        cur = con.execute("""\ 
     578            SELECT synonyms.synonym 
     579            FROM synonyms NATURAL LEFT JOIN genes 
     580            WHERE genes.gene_name=? 
     581            """, (id,)) 
     582        return map(itemgetter(0), cur) 
     583         
     584    def all_edges(self, taxid=None): 
     585        """ Return a list of all edges. 
     586        """ 
     587        con = self._db(self.taxid) 
     588        cur = con.execute(""" 
     589            SELECT links.gene_a, links.gene_b, links.weight 
     590            FROM links""") 
     591        id_to_name = self._gene_id_to_name() 
     592        return [(id_to_name[r[0]], id_to_name[r[1]], r[2]) \ 
     593                for r in cur] 
     594         
     595    def all_edges_annotated(self, taxid=None): 
     596        """ Return a list of all edges with all available annotations 
     597        """ 
     598        con = self._db(self.taxid) 
     599        cur = con.execute("""\ 
     600            SELECT links.gene_a, links.gene_b, links.weight, links.network_id 
     601            FROM links""") 
     602        gene_to_name = self._gene_id_to_name() 
     603        network_to_description = self._network_id_to_description() 
     604        res = [] 
     605        for gene_a, gene_b, w, n_id in cur: 
     606            n_desc = network_to_description[n_id] 
     607             
     608            res.append(GeneManiaInteraction(gene_to_name[gene_a], 
     609                            gene_to_name[gene_b], w, *n_desc)) 
     610        return res 
     611         
     612    def edges(self, id1): 
     613        """ Return all edges for primary id `id1`. 
     614        """         
     615        con = self._db(self.taxid) 
     616        cur = con.execute("""\ 
     617            SELECT genes1.gene_name, genes2.gene_name, links.weight 
     618            FROM genes AS genes1   
     619                JOIN links 
     620                    ON genes1.internal_id=links.gene_a 
     621                JOIN genes AS genes2 
     622                    ON genes2.internal_id=links.gene_b 
     623            WHERE genes1.gene_name=? 
     624            """, (id1,)) 
     625        res = cur.fetchall() 
     626        cur = con.execute("""\ 
     627            SELECT genes1.gene_name, genes2.gene_name, links.weight 
     628            FROM genes AS genes1   
     629                JOIN  links 
     630                    ON genes1.internal_id=links.gene_a 
     631                JOIN genes AS genes2 
     632                    ON genes2.internal_id=links.gene_b 
     633            WHERE genes2.gene_name=? 
     634            """, (id1,)) 
     635        res += cur.fetchall() 
     636         
     637        return res 
     638     
     639    def edges_annotated(self, id=None): 
     640        """ Return a list of annotated edges for primary `id`  
     641        """ 
     642        con = self._db(self.taxid) 
     643        cur = con.execute("""\ 
     644            SELECT genes1.gene_name, genes2.gene_name, links.weight, 
     645                   networks.network_name, networks.network_group, 
     646                   networks.source, networks.pubmed_id 
     647            FROM genes AS genes1 
     648                JOIN  links 
     649                    ON genes1.internal_id=links.gene_a 
     650                JOIN genes AS genes2 
     651                    ON genes2.internal_id=links.gene_b 
     652                NATURAL JOIN networks 
     653            WHERE genes1.gene_name=? 
     654            """, (id,)) 
     655        res = cur.fetchall() 
     656        cur = con.execute("""\ 
     657            SELECT genes1.gene_name, genes2.gene_name, links.weight, 
     658                   networks.network_name, networks.network_group, 
     659                   networks.source, networks.pubmed_id 
     660            FROM genes AS genes1 
     661                JOIN links 
     662                    ON genes1.internal_id=links.gene_a 
     663                JOIN genes AS genes2 
     664                    ON genes2.internal_id=links.gene_b 
     665                NATURAL JOIN networks 
     666            WHERE genes2.gene_name=? 
     667            """, (id,)) 
     668        res += cur.fetchall() 
     669        return [GeneManiaInteraction(*r) for r in res] 
     670     
     671    def search_id(self, name, taxid=None): 
     672        """ Search the database for gene name. Return a list of matching  
     673        primary ids. Use `taxid` to limit the results to a single organism. 
     674         
     675        """ 
     676        con = self._db(self.taxid) 
     677        cur = con.execute("""\ 
     678            SELECT genes.gene_name 
     679            FROM genes NATURAL JOIN synonyms 
     680            WHERE synonyms.synonym=?  
     681            """, (name,)) 
     682        return map(itemgetter(0), cur) 
     683         
     684    def _db(self, taxid=None): 
     685        """ Return an open sqlite3.Connection object.   
     686        """ 
     687        taxid = taxid or self.taxid 
     688        filename = orngServerFiles.localpath_download("PPI", 
     689                            self.SERVER_FILE.format(taxid=taxid)) 
     690        if not os.path.exists(filename): 
     691            raise ValueError("Database is missing.") 
     692         
     693        return sqlite3.connect(filename) 
     694     
     695    @lru_cache(maxsize=1) 
     696    def _gene_id_to_name(self): 
     697        """ Return a dictionary mapping internal gene ids to  
     698        primary gene identifiers. 
     699         
     700        """ 
     701        con = self._db(self.taxid) 
     702        cur = con.execute("SELECT * FROM genes") 
     703        return dict(cur) 
     704     
     705    @lru_cache(maxsize=1) 
     706    def _network_id_to_description(self): 
     707        """ Return a dictionary mapping internal network ids 
     708        to (name, group, source, pubmed id). 
     709          
     710        """ 
     711        con = self._db(self.taxid) 
     712        cur = con.execute("SELECT * FROM networks") 
     713        return dict((t[0], t[1:]) for t in cur) 
     714     
     715    ##################################### 
     716    # Data download and DB initialization 
     717    ##################################### 
     718      
     719    @classmethod 
     720    def download_data(cls, taxid=None, progress_callback=None): 
     721        """ Download the data for ``taxid`` from the GeneMANIA 
     722        website and initialize the local database. 
     723         
     724        """ 
     725        import tarfile 
     726         
     727        baseurl = "http://genemania.org/data/current/" 
     728        directory = orngServerFiles.localpath("PPI") 
     729        if taxid is None: 
     730            taxid = cls.common_taxids() 
     731         
     732        if isinstance(taxid, (list, tuple)): 
     733            taxids = taxid 
     734        else: 
     735            taxids = [taxid] 
     736        for taxid in taxids: 
     737            name = obiTaxonomy.name(taxid) 
     738            name = name.replace(" ", "_") 
     739             
     740            if progress_callback is None: 
     741                progress = True #orngServerFiles.ConsoleProgressBar("Downloading %r." % filename) 
     742            else: 
     743                progress = progress_callback 
     744             
     745            filename = name + ".tgz" 
     746            url = baseurl + "networks/" + filename     
     747            wget(url, directory=directory, progress=progress) 
     748             
     749            tgz_filename = os.path.join(directory, filename)     
     750            tgz = tarfile.open(tgz_filename) 
     751            tgz.extractall(directory) 
     752             
     753            filename = name + ".COMBINED.tgz" 
     754            url = baseurl + "precombined/" + filename 
     755            wget(url, directory=directory, progress=progress) 
     756             
     757            tgz_filename = os.path.join(directory, filename) 
     758            tgz = tarfile.open(tgz_filename) 
     759            tgz.extractall(directory) 
     760         
     761            cls.init_db([taxid]) 
     762         
     763    @classmethod 
     764    def init_db(cls, taxid=None): 
     765        """ Init the local data base. 
     766        """ 
     767        from functools import partial 
     768        directory = orngServerFiles.localpath("PPI") 
     769        pjoin = partial(os.path.join, directory) 
     770        if taxid is None: 
     771            taxid = cls.common_taxids() 
     772             
     773        if isinstance(taxid, (list, tuple)): 
     774            for tid in taxid: 
     775                cls.init_db(tid) 
     776            return 
     777                 
     778        if not isinstance(taxid, basestring): 
     779            raise ValueError("wrong taxid") 
     780             
     781#        taxid = taxids 
     782        name = obiTaxonomy.name(taxid).replace(" ", "_") 
     783        networks = csv.reader(open(pjoin(name, "networks.txt")), delimiter="\t") 
     784        networks.next() # Header 
     785        networks = list(networks) 
     786         
     787        database = pjoin(cls.SERVER_FILE.format(taxid=taxid)) 
     788        with sqlite3.connect(database) as con: 
     789            con.execute("""DROP TABLE IF EXISTS genes""") 
     790            con.execute("""DROP TABLE IF EXISTS synonyms""") 
     791            con.execute("""DROP TABLE IF EXISTS source""") 
     792            con.execute("""DROP TABLE IF EXISTS links""") 
     793            con.execute("""DROP TABLE IF EXISTS networks""") 
     794             
     795            con.execute("""DROP INDEX IF EXISTS genes_index""") 
     796            con.execute("""DROP INDEX IF EXISTS links_index_a""") 
     797            con.execute("""DROP INDEX IF EXISTS links_index_b""") 
     798             
     799            con.execute("""\ 
     800                CREATE TABLE networks  
     801                    (network_id INTEGER PRIMARY KEY ASC AUTOINCREMENT, 
     802                     network_name TEXT, 
     803                     network_group TEXT, 
     804                     source TEXT, 
     805                     pubmed_id TEXT 
     806                    )""") 
     807             
     808            con.executemany("""\ 
     809                INSERT INTO networks 
     810                VALUES (?, ?, ?, ?, ?)""", [(i, r[2], r[1], r[3], r[4]) \ 
     811                                        for i, r in enumerate(networks)]) 
     812             
     813            con.execute("""\ 
     814                CREATE TABLE genes 
     815                    (internal_id INTEGER PRIMARY KEY ASC AUTOINCREMENT, 
     816                     gene_name TEXT 
     817                    )""") 
     818             
     819            identifiers = csv.reader(open(pjoin(name, "identifier_mappings.txt"), "rb"), 
     820                                    delimiter="\t") 
     821            identifiers.next() # skip header 
     822            identifiers = list(identifiers) 
     823            genes = sorted(set(r[0] for r in identifiers)) 
     824            sources = sorted(set(r[2] for r in identifiers)) 
     825             
     826            con.executemany("""\ 
     827                INSERT INTO genes 
     828                VALUES (?, ?)""", enumerate(genes)) 
     829             
     830            con.execute("""\ 
     831            CREATE TABLE source 
     832                (source_id INTEGER PRIMARY KEY ASC AUTOINCREMENT, 
     833                 source_name TEXT 
     834                )""") 
     835             
     836            con.executemany("""\ 
     837                INSERT INTO source 
     838                VALUES (?, ?)""", enumerate(sources)) 
     839             
     840            con.execute("""\ 
     841                CREATE TABLE synonyms 
     842                    (internal_id INTEGER REFERENCES genes (internal_id), 
     843                     synonym TEXT, 
     844                     source_id INT REFERENCES source (source_id) 
     845                    )""") 
     846             
     847            gene_to_id = dict((g, i) for i, g in enumerate(genes)) 
     848            source_to_id = dict((s, i) for i, s in enumerate(sources)) 
     849            con.executemany("""\ 
     850                INSERT INTO synonyms 
     851                VALUES (?, ?, ?)""", [(gene_to_id[r[0]], r[1], source_to_id[r[2]])\ 
     852                                       for r in identifiers]) 
     853             
     854            con.execute("""\ 
     855                CREATE TABLE links 
     856                    (gene_a INTEGER REFERENCES genes (internal_id), 
     857                     gene_b INTEGER REFERENCES genes (internal_id), 
     858                     network_id INTEGER REFERENCES networks (network_id), 
     859                     weight REAL 
     860                     -- , PRIMARY KEY (gene_a, gene_b, network_id) 
     861                    )""") 
     862             
     863            for i, (filename, group, _, _, _) in enumerate(networks): 
     864                nf  = open(pjoin(name, filename), "rb") 
     865                interactions = csv.reader(nf, delimiter="\t") 
     866                interactions.next() # skip header 
     867                con.executemany("""\ 
     868                    INSERT INTO links  
     869                    VALUES (?, ?, ?, ?)""", 
     870                    [(gene_to_id[r[0]], gene_to_id[r[1]], i, float(r[2])) \ 
     871                     for r in interactions] 
     872                ) 
     873                 
     874            # Add special combined network entry 
     875            combined_id = len(networks) 
     876            con.execute("""\ 
     877                INSERT INTO networks 
     878                VALUES (?, ?, ?, ?, ?)""",  
     879                (combined_id, "BP_COMBINING", "COMBINED", "GeneMANIA", "")) 
     880             
     881            # Add the combined network links. 
     882            combined = open(pjoin(name + ".COMBINED", "COMBINED.DEFAULT_NETWORKS.BP_COMBINING.txt"), "rb") 
     883            combined = csv.reader(combined, delimiter="\t") 
     884            combined.next() 
     885            con.executemany("""\ 
     886                INSERT INTO links 
     887                VALUES (?, ?, ?, ?)""", 
     888                    ((gene_to_id[r[0]], gene_to_id[r[1]], combined_id, float(r[2])) \ 
     889                     for r in combined)) 
     890             
     891             
     892            con.execute(""" 
     893                CREATE VIEW IF NOT EXISTS links_annotated 
     894                AS SELECT genes1.gene_name AS gene_name_a, 
     895                          genes2.gene_name AS gene_name_b, 
     896                          links.weight, 
     897                          networks.network_name, 
     898                          networks.network_group, 
     899                          networks.source, 
     900                          networks.pubmed_id 
     901                   FROM  genes AS genes1 
     902                        JOIN links  
     903                              ON genes1.internal_id=links.gene_a 
     904                        JOIN genes AS genes2 
     905                              ON links.gene_b=genes2.internal_id 
     906                        JOIN networks 
     907                              ON links.network_id=networks.network_id 
     908                    """) 
     909             
     910             
     911            con.execute("""\ 
     912                CREATE INDEX IF NOT EXISTS genes_index ON genes (gene_name) 
     913                """) 
     914            con.execute("""\ 
     915                CREATE INDEX IF NOT EXISTS links_index_a ON links (gene_a) 
     916                """) 
     917            con.execute("""\ 
     918                 CREATE INDEX IF NOT EXISTS links_index_b ON links (gene_b) 
     919                """) 
     920         
     921             
    391922if __name__ == "__main__": 
    392923    retrieve("9606", [ 'MRE11A', 'RAD51', 'MLH1', 'MSH2', 'DMC1', 'RAD51AP1', 'RAD50', 'MSH6', 'XRCC3', 'PCNA', 'XRCC2' ]) 
Note: See TracChangeset for help on using the changeset viewer.