Pandas API for Gene Set Enrichment Analysis in Python (GSEApy, cudaGSEA, GSEA)
- aims to provide a unified API for various GSEA implementations; uses pandas DataFrames and a hierarchy of Pythonic classes.
- file exports (exporting input for GSEA) use low-level numpy functions and are much faster than in pandas
- aims to allow researchers to easily compare different implementations of GSEA, and to integrate those in projects which require high-performance GSEA (e.g. massive screening for drug-repositioning)
- provides useful utilities for work with GMT files, or gene sets and pathways in general in Python
To install the API use:
pip3 install gsea_apiSee below for the instructions on installation of specific GSEA implementations.
frompandasimportread_tablefromgsea_api.expression_setimportExpressionSetfromgsea_api.gseaimportGSEADesktopfromgsea_api.molecular_signatures_dbimportGeneSetsreactome_pathways=GeneSets.from_gmt('ReactomePathways.gmt')
gsea=GSEADesktop()
design= ['Disease', 'Disease', 'Disease', 'Control', 'Control', 'Control']
matrix=read_table('expression_data.tsv', index_col='Gene')
result=gsea.run(
# note: contrast() is not necessary in this simple caseExpressionSet(matrix, design).contrast('Disease', 'Control'),
reactome_pathways,
metric='Signal2Noise',
permutations=1000
)Where expression_data.tsv is in the following format:
Gene Patient_1 Patient_2 Patient_3 Patient_4 Patient_5 Patient_6
TACC2 0.2 0.1 0.4 0.6 0.7 2.1
TP53 2.3 0.2 2.1 2.0 0.3 0.6
Molecular Signatures Database (MSigDB) can be downloaded from the Broad Institute GSEA website. It provides expert-curated gene set collections, as well as curated subset of pathway databases (Reactome, KEGG, Biocarta, Gene Ontology) trimmed to remove redundant, overlapping and and otherwise little-value terms (if needed).
You can download all the pathways collections at once (search for ZIPped MSigDB on the download page). After downloading and un-zipping (e.g., to a local directory named msigdb), you can access the gene sets from MSigDB with:
fromgsea_api.molecular_signatures_dbimportMolecularSignaturesDatabasemsigdb=MolecularSignaturesDatabase('msigdb', version=7.1)
msigdb.gene_setsmsigdb.gene_sets returns a list of dictionaries describing auto-detected pathways:
[
{'name': 'c1.all', 'id_type': 'symbols'},
{'name': 'c1.all', 'id_type': 'entrez'},
{'name': 'c2.cp.reactome', 'id_type': 'symbols'},
{'name': 'c2.cp.reactome', 'id_type': 'entrez'}
# etc..
]Information about the location on disk and version are available in msigdb.path and msigdb.version.
msigdb.load loads the specific collection into a GeneSets object:
>kegg_pathways=msigdb.load('c2.cp.kegg', 'symbols')
>print(kegg_pathways)
<GeneSets'c2.cp.kegg'with186genesets>This object can be passed to any of the supported GSEA implementations; please see below for a detailed description of the GeneSets object.
GeneSets represents a collection of sets of genes, where each set is represented as GeneSet object.
You can check the number of sets contained within a collection with:
>len(kegg_pathways)
186The gene sets are accessible with gene_sets (tuple) and gene_sets_by_name (dict) properties:
>kegg_pathways.gene_sets[:2]
(<GeneSet'KEGG_TIGHT_JUNCTION'with132genes>, <GeneSet'KEGG_RNA_DEGRADATION'with59genes>)
>kegg_pathways.gene_sets_by_name
{
'KEGG_TIGHT_JUNCTION': <GeneSet'KEGG_TIGHT_JUNCTION'with132genes>,
'KEGG_RNA_DEGRADATION': <GeneSet'KEGG_RNA_DEGRADATION'with59genes># etc.
}Sometimes only a subset of genes is measured in an experiment. You can remove gene sets which do not contain any of the measured genes from the collection:
>measured_genes= {'APOE', 'CYB5R1', 'FCER1G', 'PVR', 'HK2'}
>measured_subset=kegg_pathways.subset(measured_genes)
>print(measured_subset)
<GeneSetswith12genesets>The skipped gene sets are accessible in measured_subset.empty_gene_sets for inspection.
>kegg_pathways.trim(min_genes=10, max_genes=20)
<GeneSetswith21genesets>defprettify_kegg_name(gene_set):
returngene_set.name.replace('KEGG_', '').replace('_', ' ')
kegg_pathways_pretty=kegg_pathways.format_names(prettify_kegg_name)
kegg_pathways_pretty.gene_sets[:2]
# (<GeneSet 'TIGHT JUNCTION' with 132 genes>, <GeneSet 'RNA DEGRADATION' with 59 genes>)For MSigDB 7.4+:
defpretty_reactome_name(gene_set):
returngene_set.metadata['DESCRIPTION_BRIEF']
reactome_pathways_pretty=reactome_pathways.format_names(pretty_reactome_name)
reactome_pathways_pretty.gene_sets[:2]
#Other properties and methods offered by GeneSets include:
all_genes: return a set of all genes which are covered by the gene sets in the collectionname: the name of the collectionto_frame()return a pandasDataFramedescribing membership of the genes (gene sets = rows, genes = columns), which can be used for UpSet visualisation (e.g. with ComplexUpset)to_gmt(path: str)exports the gene set to a GMT (Gene Matrix Transposed) file
Following GSEA implementations are supported:
Login/register on the official GSEA website and download the gsea_3.0.jar file (or a newer version).
Provide the location of the downloaded file to GSEADesktop() using gsea_jar_path argument, e.g.:
gsea=GSEADesktop(gsea_jar_path='downloads/gsea_3.0.jar')To use gsea.py please install it with:
pip3 install gseapy
Use it with:
fromgsea_api.gseaimportGSEApygsea=GSEApy()Please clone this fork of cudaGSEA and compile the binary version:
git clone https://github.com/krassowski/cudaGSEA
cd cudaGSEA/cudaGSEA/src/
# if on Ubuntu:# sudo apt install nvidia-cuda-toolkit# whereis nvccexport CUDA_HOME=/usr
export R_INC=/usr/share/R/include
export RCPP_INC=/usr/local/lib/R/site-library/Rcpp/include
make cudaGSEAdepending on your GPU and drivers you may see Unsupported gpu architecture 'compute_20' error; simply edit Makefile removing -gencode arch=compute_20,code=compute_20 (see this askUbuntu post)
You can also try to use the original version, which does not implement FDR calculations.
Use it with:
fromgsea_api.gseaimportcudaGSEA# CPU implementation can be used with use_cpu=Truegsea=cudaGSEA(fdr='full', use_cpu=False, path='cudaGSEA/cudaGSEA/src/cudaGSEA')Please also cite the authors of the wrapped tools that you use.
The initial version of this code was written for a Master thesis project at Imperial College London.