pheno_prerank_enrich() Enrichment: Gene Interaction List with Rank Scores for Phenotypes
Step 1: Prepare rank_df
Begin by preparing a list of gene interactions along with their preranked scores. Store this information in a DataFrame called rank_df.
- The index of
rank_dfshould be the gene interaction symbols. e.g.Oas1a_Ifit1 - The columns of
rank_dfshould correspond to one target phenotype (e.g., GK vs.).
This format allows pheno_prerank_enrich to perform enrichment analysis using phenotype-level ranked interaction-level statistics.
In this tutorial, we use gene interaction data from the GSE13268 dataset, which includes two phenotypes:
-
GK (Goto-Kakizaki): An inbred rat model commonly used for diabetes research. GK rats exhibit a polygenic form of diabetes, closely mirroring human disease characteristics such as hyperglycemia, impaired glucose tolerance, and insulin resistance. This makes them a valuable model for studying human type 2 diabetes.
-
WKY (Wistar-Kyoto): A standard inbred laboratory rat strain often used as a control. In this study, WKY rats serve as the reference group for comparison with GK rats, enabling the investigation of insulin resistance and diabetes progression.
For each rat sample, gene–gene interactions and their corresponding entropy-based Critical Transition (CT) scores were precomputed by NIEE. A differential analysis was then performed between the GK and WKY groups. From this, we identified gene interactions with higher fluctuation in GK rats, ranked by their z-scores from the differential test (grea/data/GSE13268_GK.csv).
%load_ext autoreload
import pandas as pd
rank_df = pd.read_csv("grea/data/GSE13268_GK.csv", index_col=0)
rank_df.head()
Step 2: Preparing Gene Set Libraries
There are several ways to prepare gene set libraries for use in GREA:
Option 1: Use Built-in Libraries
Simply specify the libraries you're interested in as a list. For example:
libraries = ['KEGG_2019_Mouse', 'WikiPathways_2024_Mouse']
You can use the grea.library.list_libraries() function to view all available pathway libraries included in GREA.
Option 2: Load from GMT File
You can load external gene set libraries from .gmt files using:
libraries = read_gmt('your_library_file.gmt')
Option 3: Define a Custom Library
Create your own gene set library using a Python dictionary, where each key is a pathway name and the corresponding value is a list of genes:
libraries = {
'term1': ['Oas1a', 'Ifit1'],
'term2': ['Gstm4', 'Oas1a']
}
from grea.library import list_libraries
print(list_libraries())
Step 3: Run Enrichment
To perform enrichment analysis, call the grea.pheno_prerank_enrich(rank_df, libraries) function. You can customize the analysis using the following arguments:
prob_method: Method for p-value calculation. Currently supports'perm'for permutation-based testing.n_perm: Number of permutations to use for estimating the null distribution.sig_sep: The delimiter used to separate gene names in an interaction string (e.g., setsig_sep='_'for interactions likeOas1a_Ifit1).
The function returns a GREA object containing all enrichment results, including enrichment scores and statistical significance for each library term.
%autoreload
from grea import grea
libraries = ['WikiPathways_2024_Mouse', 'Mouse_Gene_Atlas']
n_perm = 10
prob_method = 'perm'
sig_sep = '_'
obj = grea.pheno_prerank_enrich(rank_df, libraries,n_perm=n_perm, prob_method=prob_method, sig_sep=sig_sep)
obj
Step 4: Check Enrichment Results
The GREA object stores all enrichment results, including enrichment scores and statistical significance for each library term. GREA supports three types of enrichment scores, each reflecting a different scoring strategy:
'KS-ES': Kolmogorov–Smirnov-based Enrichment Score, capturing the peak deviation between hit and miss distributions.'KS-ESD': KS-based enrichment Score Difference, the sum of the maximum positive and negative deviations from the running score.'RC-AUC': Area Under the Recovery Curve, summarizing early enrichment of target genes along the ranking.
You can select the appropriate metric depending on your analysis goal or data characteristics.
To retrieve the enrichment results as a long DataFrame, use the get_enrich_results(metric) function.
%autoreload
df = obj.get_enrich_results(metric='KS-ES')
df.head()
%autoreload
df = obj.get_enrich_results(metric='KS-ESD')
df.head()
%autoreload
df = obj.get_enrich_results(metric='RC-AUC')
df.head()
Step 5: Visualize Enrichment Results
To visualize the enrichment results, use the pl_running_sum(metric, term, pheno_id) function by specifying the desired metric, term, and target phenotype.
%autoreload
term = 'Mouse_Gene_Atlas|uterus'
pheno_id = 'GK vs. WKY'
fig = obj.pl_running_sum('KS-ES', term, pheno_id)
fig
%autoreload
term = 'Mouse_Gene_Atlas|uterus'
pheno_id = 'GK vs. WKY'
fig = obj.pl_running_sum('KS-ESD', term, pheno_id)
fig
%autoreload
term = 'Mouse_Gene_Atlas|uterus'
pheno_id = 'GK vs. WKY'
fig = obj.pl_running_sum('RC-AUC', term, pheno_id)
fig
