Gene set selections in Python

Overview

The gesel package provides a Python interface to the Gesel database for client-side gene set searches. The idea is to use HTTP range requests to serve the gene-to-set mappings to a client without downloading the entire database or implementing custom server logic. In this manner, we can execute a variety of interesting gene set queries with high scalability across users and minimal backend maintenance.

Quick start

To get started, install the package from PyPI:

pip install gesel

Then we can find overlaps between our genes of interest and the gene sets in the Gesel database. (Note that the exact numbers are subject to change, pending updates to the version of the Gesel database.)

import gesel
my_genes = ["SNAP25", "NEUROD6", "GAD1", "GAD2", "RELN"]

# First, mapping our gene names to Gesel's internal gene indices.
gene_idx = gesel.search_genes("9606", my_genes) # list of lists of gene indices.
print(gene_idx)
## [[4639], [12767], [1758], [1759], [3912]]
gene_idx = sum(gene_idx, []) # collapsing it to a list of integers, for simplicity.
print(gesel.fetch_all_genes("9606")[gene_idx,:]) # double-checking that we got it right.
## BiocFrame with 5 rows and 3 columns
##          symbol    entrez             ensembl
##          <list>    <list>              <list>
## [0]  ['SNAP25']  ['6616'] ['ENSG00000132639']
## [1] ['NEUROD6'] ['63974'] ['ENSG00000164600']
## [2]    ['GAD1']  ['2571'] ['ENSG00000128683']
## [3]    ['GAD2']  ['2572'] ['ENSG00000136750']
## [4]    ['RELN']  ['5649'] ['ENSG00000189056']

# Now finding all sets with one or more overlaps to `my_genes`.
overlaps, present = gesel.find_overlapping_sets("9606", gene_idx, counts_only=False)
print(overlaps) # set index and the identities of overlapping genes.
## BiocFrame with 1163 rows and 2 columns
##           set                     genes
##        <list>                    <list>
##    [0]   2420  [4639, 1758, 1759, 3912]
##    [1]   2521  [4639, 1758, 1759, 3912]
##    [2]  21748 [4639, 12767, 1758, 1759]
##           ...                       ...
## [1160]  40562                    [3912]
## [1161]  40597                    [3912]
## [1162]  40599                    [3912]

# Actually getting the identities of the top sets:
set_info = gesel.fetch_some_sets("9606", overlaps["set"][:10])
print(set.info)
## BiocFrame with 10 rows and 5 columns
##                         name             description   size collection number
##                       <list>                  <list> <list>     <list> <list>
## [0]               GO:0005737               cytoplasm   5010          0   2420
## [1]               GO:0005886         plasma membrane   4778          0   2521
## [2]  BLALOCK_ALZHEIMERS_D... http://www.gsea-msig...   1248          2   2515
## [3]  MANNO_MIDBRAIN_NEURO... http://www.gsea-msig...   1106         14     92
## [4]               GO:0005515         protein binding  12505          0   2310
## [5]               GO:0007268 chemical synaptic tr...    248          0   3331
## [6]  MIKKELSEN_MEF_HCP_WI... http://www.gsea-msig...    591          2   1896
## [7]  REACTOME_NEUROTRANSM... http://www.gsea-msig...     51          4     29
## [8]  REACTOME_TRANSMISSIO... http://www.gsea-msig...    270          4     32
## [9] REACTOME_NEURONAL_SYSTEM http://www.gsea-msig...    411          4     33

# As well as the collections from which they were derived.
print(gesel.fetch_some_collections("9606", [0, 2, 4, 14]))
## BiocFrame with 4 rows and 6 columns
##                       title             description maintainer                  source  start   size
##                      <list>                  <list>     <list>                  <list> <list> <list>
## [0]           Gene ontology Gene sets defined fr...  Aaron Lun https://github.com/L...      0  18933
## [1] MSigDB chemical and ... Gene sets that repre...  Aaron Lun https://github.com/L...  19233   3405
## [2] MSigDB canonical pat... Reactome gene sets a...  Aaron Lun https://github.com/L...  22834   1654
## [3] MSigDB cell type sig... Gene sets that conta...  Aaron Lun https://github.com/L...  39774    830

Check out the reference documentation for more details.

Searching on text

We can also search for gene sets based on the text in their names or descriptions.

chits = gesel.search_set_text("9606", "cancer")
print(gesel.fetch_some_sets("9606", chits[:10]))
## BiocFrame with 10 rows and 5 columns
##                        name             description   size collection number
##                      <list>                  <list> <list>     <list> <list>
## [0] SOGA_COLORECTAL_CANC... http://www.gsea-msig...     71          2      1
## [1] SOGA_COLORECTAL_CANC... http://www.gsea-msig...     82          2      2
## [2] WATANABE_RECTAL_CANC... http://www.gsea-msig...    113          2     64
## [3]  LIU_PROSTATE_CANCER_UP http://www.gsea-msig...     99          2     66
## [4] BERTUCCI_MEDULLARY_V... http://www.gsea-msig...    207          2     68
## [5] WATANABE_COLON_CANCE... http://www.gsea-msig...     29          2     78
## [6] WATANABE_COLON_CANCE... http://www.gsea-msig...     69          2     79
## [7] SOTIRIOU_BREAST_CANC... http://www.gsea-msig...     53          2     81
## [8] CHARAFE_BREAST_CANCE... http://www.gsea-msig...     52          2    124
## [9] DOANE_BREAST_CANCER_... http://www.gsea-msig...     33          2    125

ihits = gesel.search_set_text("9606", "innate immun*")
print(gesel.fetch_some_sets("9606", ihits[:10]))
## BiocFrame with 10 rows and 5 columns
##                        name             description   size collection number
##                      <list>                  <list> <list>     <list> <list>
## [0] REACTOME_INNATE_IMMU... http://www.gsea-msig...   1118          4    229
## [1] REACTOME_REGULATION_... http://www.gsea-msig...     15          4    552
## [2] REACTOME_SARS_COV_1_... http://www.gsea-msig...     41          4   1562
## [3] REACTOME_SARS_COV_2_... http://www.gsea-msig...    126          4   1580
## [4] WP_SARSCOV2_B117_VAR... http://www.gsea-msig...      9          5    185
## [5] WP_SARS_CORONAVIRUS_... http://www.gsea-msig...     31          5    189
## [6] WP_SARSCOV2_INNATE_I... http://www.gsea-msig...     66          5    318
## [7] WP_PATHWAYS_OF_NUCLE... http://www.gsea-msig...     16          5    729
## [8]              GO:0002218 activation of innate...     32          0    813
## [9]              GO:0002220 innate immune respon...      1          0    814

thits = gesel.search_set_text("9606", "cd? t cell")
print(gesel.fetch_some_sets("9606", thits[:10]))
##                        name             description   size collection number
##                      <list>                  <list> <list>     <list> <list>
## [0] HOFT_CD4_POSITIVE_AL... http://www.gsea-msig...     47         13     49
## [1] HOFT_CD4_POSITIVE_AL... http://www.gsea-msig...     28         13    131
## [2] HOFT_CD4_POSITIVE_AL... http://www.gsea-msig...     40         13    204
## [3] HOFT_CD4_POSITIVE_AL... http://www.gsea-msig...     16         13    252
## [4] HOFT_CD4_POSITIVE_AL... http://www.gsea-msig...     24         13    296
## [5] HOFT_CD4_POSITIVE_AL... http://www.gsea-msig...     42         13    316
## [6] HOFT_CD4_POSITIVE_AL... http://www.gsea-msig...     41         13    320
## [7] HOFT_CD4_POSITIVE_AL... http://www.gsea-msig...      6         13    323
## [8] QI_CD4_POSITIVE_ALPH... http://www.gsea-msig...      9         13    327
## [9] QI_NAIVE_T_CELL_ZOST... http://www.gsea-msig...      7         13    328

Users can construct powerful queries by intersecting the sets recovered from search_set_text() with those from find_overlapping_sets().

import biocutils
cancer_sets = biocutils.intersect(chits, overlaps["set"])
info = gesel.fetch_some_sets("9606", cancer_sets)
m = biocutils.match(cancer_sets, overlaps["set"])
info = info.set_column("count", [len(overlaps["genes"][r]) for r in m])

# We'll just use the proportion of enriched genes for ranking here;
# a more sophisticated analysis might compute a hypergeometric p-value.
prop = [info["count"][i] / info["size"][i] for i in range(info.shape[0])]
ordered = info[biocutils.order(prop, decreasing=True),:]
print(ordered)
## BiocFrame with 14 rows and 6 columns
##                         name             description   size collection number  count
##                       <list>                  <list> <list>     <list> <list> <list>
##  [0] LOPES_METHYLATED_IN_... http://www.gsea-msig...     28          2    903      1
##  [1] SCHLESINGER_H3K27ME3... http://www.gsea-msig...     28          2   2864      1
##  [2] WATANABE_COLON_CANCE... http://www.gsea-msig...     29          2     78      1
##                          ...                     ...    ...        ...    ...    ...
## [11] ACEVEDO_LIVER_CANCER_DN http://www.gsea-msig...    540          2   3150      1
## [12] SMID_BREAST_CANCER_L... http://www.gsea-msig...    587          2    991      1
## [13] LIU_OVARIAN_CANCER_T... http://www.gsea-msig...   1713          2   2003      1

Fetching all data

gesel is designed around partial extraction from the database files, but it may be more efficient to pull all of the data into memory at once. This is most useful for the gene set details, which can be retrieved en masse:

set_info = gesel.fetch_all_sets("9606")
print(set_info)
## BiocFrame with 40654 rows and 5 columns
##                            name             description   size collection number
##                          <list>                  <list> <list>     <list> <list>
##     [0]              GO:0000002 mitochondrial genome...     11          0      0
##     [1]              GO:0000003            reproduction      4          0      1
##     [2]              GO:0000009 alpha-1,6-mannosyltr...      2          0      2
##                             ...                     ...    ...        ...    ...
## [40651] HALLMARK_KRAS_SIGNAL... http://www.gsea-msig...    200         15     47
## [40652] HALLMARK_KRAS_SIGNAL... http://www.gsea-msig...    200         15     48
## [40653] HALLMARK_PANCREAS_BE... http://www.gsea-msig...     40         15     49

The set indices returned by other functions like find_overlapping_sets() can then be used to directly subset the set_info data frame by row. In fact, calling fetch_some_sets() after fetch_all_sets() will automatically use the data frame created by the latter, instead of attempting another request to the database. The same approach can be used to extract collection information, via fetch_all_collections(); gene set membership, via fetch_genes_for_all_sets(); and the sets containing each gene, via fetch_sets_for_all_genes().

Advanced use

gesel uses a lot of in-memory caching to reduce the number of requests to the database files within a single Python session. On rare occasions, the cache may become outdated, e.g., if the database files are updated while an Python session is running. Users can prompt gesel to re-acquire all data by flusing the cache:

gesel.flush_memory_cache()

Applications can specify their own functions for obtaining files (or byte ranges thereof) by passing a custom config= in each gesel function. For example, on a shared HPC filesystem, we could point gesel towards a directory of Gesel database files. This provides a performant alternative to HTTP requests for an institutional collection of gene sets.