<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 GSEAPY Example
Examples to use GSEApy inside python console
Step1:  Check gseapy version
Step2:  1.  (Optional) Convert IDs Using Biomart API
Don't use this if you don't know Biomart
python
    &gt;&gt;&gt; from gseapy.parser import Biomart 
    &gt;&gt;&gt; bm = Biomart()
    &gt;&gt;&gt; ## view validated marts
    &gt;&gt;&gt; marts = bm.get_marts()
    &gt;&gt;&gt; ## view validated dataset
    &gt;&gt;&gt; datasets = bm.get_datasets(mart='ENSEMBL_MART_ENSEMBL')
    &gt;&gt;&gt; ## view validated attributes
    &gt;&gt;&gt; attrs = bm.get_attributes(dataset='hsapiens_gene_ensembl') 
    &gt;&gt;&gt; ## view validated filters
    &gt;&gt;&gt; filters = bm.get_filters(dataset='hsapiens_gene_ensembl')
    &gt;&gt;&gt; ## query results
    &gt;&gt;&gt; queries = ['ENSG00000125285','ENSG00000182968'] # need to be a python list
    &gt;&gt;&gt; results = bm.query(dataset='hsapiens_gene_ensembl', 
                           attributes=['ensembl_gene_id', 'external_gene_name', 'entrezgene_id', 'go_id'],
                           filters={'ensemble_gene_id'
Step3:  See all supported enrichr library names  
Select database from { 'Human', 'Mouse', 'Yeast', 'Fly', 'Fish', 'Worm' }
Enrichr library could be used for gsea, ssgsea, and prerank, too
Step4:  2.1  Assign enrichr with pd.Series,  pd.DataFrame, or list object
2.1.1 gene_sets support list, str.
Multi-libraries names supported, separate each name by comma or input a list. 
For example
Step5:  2.1.2 Local mode of GO analysis
If input a .gmt file or gene_set dict object, enrichr runs local.
You have to specify the background genes, if local mode used
For example
Step6:  2.1.3 Plotting
Step7:  2.2 Command line usage
You may also want to use enrichr in command line  
the option -v will print out the progress of your job
Step8:  3. Prerank example
3.1 Assign prerank()  with a pd.DataFrame, pd.Series , or a txt file
Do not include header in your gene list !
GSEApy will skip any data after "#".
Only contains two columns, or one cloumn with gene_name indexed when assign a DataFrame to prerank
Step9:  Leading edge genes save to the final output results
Step10:  3.2 How to generate your GSEA plot inside python console
Visualize it using gseaplot  
Make sure that ofname is not None, if you want to save your figure to the disk
Step11:  3) Command line usage
You may also want to use prerank in command line
Step12:  4. GSEA Example
4.1  Assign gsea() with a pandas DataFrame, .gct format file, or a text file
and cls with a list object or just .cls format file
Step13:  4.2 Show the gsea plots
The gsea module will generate heatmap for genes in each gene sets in the backgroud.
But if you need to do it yourself, use the code below
Step14:  4.3 Command line usage
You may also want to use gsea in command line
Step15:  5. Single Sample GSEA example
Note
Step16:  5.2 Access Enrichment Score (ES) and NES
results save to two attribute
Step17:  Note
Step18:  3) command line usage of single sample gsea
Step19:  6. Replot Example
6.1 locate your directory
notes
Step20:  6.2 command line usage of replot
<ASSISTANT_TASK:>
Python Code:
# %matplotlib inline
# %config InlineBackend.figure_format='retina' # mac
# %load_ext autoreload
# %autoreload 2
import pandas as pd
import gseapy as gp
import matplotlib.pyplot as plt
Explanation: GSEAPY Example
Examples to use GSEApy inside python console
End of explanation
gp.__version__
Explanation: Check gseapy version
End of explanation
# read in an example gene list
gene_list = pd.read_csv("./tests/data/gene_list.txt",header=None, sep="\t")
gene_list.head()
# convert dataframe or series to list
glist = gene_list.squeeze().str.strip().tolist()
print(glist[:10])
Explanation: 1.  (Optional) Convert IDs Using Biomart API
Don't use this if you don't know Biomart
python
    &gt;&gt;&gt; from gseapy.parser import Biomart 
    &gt;&gt;&gt; bm = Biomart()
    &gt;&gt;&gt; ## view validated marts
    &gt;&gt;&gt; marts = bm.get_marts()
    &gt;&gt;&gt; ## view validated dataset
    &gt;&gt;&gt; datasets = bm.get_datasets(mart='ENSEMBL_MART_ENSEMBL')
    &gt;&gt;&gt; ## view validated attributes
    &gt;&gt;&gt; attrs = bm.get_attributes(dataset='hsapiens_gene_ensembl') 
    &gt;&gt;&gt; ## view validated filters
    &gt;&gt;&gt; filters = bm.get_filters(dataset='hsapiens_gene_ensembl')
    &gt;&gt;&gt; ## query results
    &gt;&gt;&gt; queries = ['ENSG00000125285','ENSG00000182968'] # need to be a python list
    &gt;&gt;&gt; results = bm.query(dataset='hsapiens_gene_ensembl', 
                           attributes=['ensembl_gene_id', 'external_gene_name', 'entrezgene_id', 'go_id'],
                           filters={'ensemble_gene_id': queries})
2. Enrichr Example
End of explanation
names = gp.get_library_name() # default: Human
names[:10]
yeast = gp.get_library_name(organism='Yeast') 
yeast[:10]
Explanation: See all supported enrichr library names  
Select database from { 'Human', 'Mouse', 'Yeast', 'Fly', 'Fish', 'Worm' }
Enrichr library could be used for gsea, ssgsea, and prerank, too
End of explanation
# run enrichr
# if you are only intrested in dataframe that enrichr returned, please set no_plot=True
# list, dataframe, series inputs are supported
enr = gp.enrichr(gene_list="./tests/data/gene_list.txt", 
                 gene_sets=['KEGG_2016','KEGG_2013'], 
                 organism='Human', # don't forget to set organism to the one you desired! e.g. Yeast
                 description='test_name',
                 outdir='test/enrichr_kegg', 
                 # no_plot=True,
                 cutoff=0.5 # test dataset, use lower value from range(0,1)
                )
# obj.results stores all results
enr.results.head(5)
Explanation: 2.1  Assign enrichr with pd.Series,  pd.DataFrame, or list object
2.1.1 gene_sets support list, str.
Multi-libraries names supported, separate each name by comma or input a list. 
For example:
python
    # gene_list
    gene_list="./data/gene_list.txt", 
    gene_list=glist
    # gene_sets
    gene_sets='KEGG_2016'  
    gene_sets='KEGG_2016,KEGG_2013'
    gene_sets=['KEGG_2016','KEGG_2013']
End of explanation
enr2 = gp.enrichr(gene_list="./tests/data/gene_list.txt", 
                 # or gene_list=glist
                 description='test_name', 
                 gene_sets="./tests/data/genes.gmt",
                 background='hsapiens_gene_ensembl', # or the number of genes, e.g 20000
                 outdir='test/enrichr_kegg2', 
                 cutoff=0.5, # only used for testing.
                 verbose=True)
enr2.results.head(5)
Explanation: 2.1.2 Local mode of GO analysis
If input a .gmt file or gene_set dict object, enrichr runs local.
You have to specify the background genes, if local mode used
For example:
python
    gene_sets="./data/genes.gmt",
    gene_sets={'A':['gene1', 'gene2',...],
               'B':['gene2', 'gene4',...],
               ...}
End of explanation
# simple plotting function
from gseapy.plot import barplot, dotplot
# to save your figure, make sure that ``ofname`` is not None
barplot(enr.res2d,title='KEGG_2013',)
# to save your figure, make sure that ``ofname`` is not None
dotplot(enr.res2d, title='KEGG_2013',cmap='viridis_r')
Explanation: 2.1.3 Plotting
End of explanation
# !gseapy enrichr -i ./data/gene_list.txt \
#                 --ds BP2017 \
#                 -g GO_Biological_Process_2017 \
#                 -v -o test/enrichr_BP
Explanation: 2.2 Command line usage
You may also want to use enrichr in command line  
the option -v will print out the progress of your job
End of explanation
rnk = pd.read_csv("./tests/data/edb/gsea_data.gsea_data.rnk", header=None, sep="\t")
rnk.head()
# run prerank
# enrichr libraries are supported by prerank module. Just provide the name
# use 4 process to acceralate the permutation speed
# note: multiprocessing may not work on windows
pre_res = gp.prerank(rnk=rnk, gene_sets='KEGG_2016', 
                     processes=4,
                     permutation_num=100, # reduce number to speed up testing
                     outdir='test/prerank_report_kegg', format='png', seed=6)
Explanation: 3. Prerank example
3.1 Assign prerank()  with a pd.DataFrame, pd.Series , or a txt file
Do not include header in your gene list !
GSEApy will skip any data after "#".
Only contains two columns, or one cloumn with gene_name indexed when assign a DataFrame to prerank
End of explanation
#access results through obj.res2d attribute or obj.results
pre_res.res2d.sort_index().head()
Explanation: Leading edge genes save to the final output results
End of explanation
# extract geneset terms in res2d
terms = pre_res.res2d.index
terms
## easy way
from gseapy.plot import gseaplot
# to save your figure, make sure that ofname is not None
gseaplot(rank_metric=pre_res.ranking, term=terms[0], **pre_res.results[terms[0]])
# save figure
# gseaplot(rank_metric=pre_res.ranking, term=terms[0], ofname='your.plot.pdf', **pre_res.results[terms[0]])
Explanation: 3.2 How to generate your GSEA plot inside python console
Visualize it using gseaplot  
Make sure that ofname is not None, if you want to save your figure to the disk
End of explanation
# ! gseapy prerank -r temp.rnk -g temp.gmt -o prerank_report_temp
Explanation: 3) Command line usage
You may also want to use prerank in command line
End of explanation
phenoA, phenoB, class_vector =  gp.parser.gsea_cls_parser("./tests/data/P53.cls")
#class_vector used to indicate group attributes for each sample
print(class_vector)
gene_exp = pd.read_csv("./tests/data/P53.txt", sep="\t")
gene_exp.head()
print("positively correlated: ", phenoA)
print("negtively correlated: ", phenoB)
# run gsea
# enrichr libraries are supported by gsea module. Just provide the name
gs_res = gp.gsea(data=gene_exp, # or data='./P53_resampling_data.txt'
                 gene_sets='KEGG_2016', # enrichr library names
                 cls= './tests/data/P53.cls', # cls=class_vector
                 # set permutation_type to phenotype if samples >=15
                 permutation_type='phenotype', 
                 permutation_num=100, # reduce number to speed up test
                 outdir=None,  # do not write output to disk
                 no_plot=True, # Skip plotting
                 method='signal_to_noise',
                 processes=4, seed= 7,
                 format='png')
#access the dataframe results throught res2d attribute
gs_res.res2d.sort_index().head()
Explanation: 4. GSEA Example
4.1  Assign gsea() with a pandas DataFrame, .gct format file, or a text file
and cls with a list object or just .cls format file
End of explanation
from gseapy.plot import gseaplot, heatmap
terms = gs_res.res2d.index
# Make sure that ``ofname`` is not None, if you want to save your figure to disk
gseaplot(gs_res.ranking, term=terms[0], **gs_res.results[terms[0]])
# plotting heatmap
genes = gs_res.res2d.genes[0].split(";")
# Make sure that ``ofname`` is not None, if you want to save your figure to disk
heatmap(df = gs_res.heatmat.loc[genes], z_score=0, title=terms[0], figsize=(18,6))
Explanation: 4.2 Show the gsea plots
The gsea module will generate heatmap for genes in each gene sets in the backgroud.
But if you need to do it yourself, use the code below
End of explanation
# !gseapy gsea -d ./data/P53_resampling_data.txt \
#              -g KEGG_2016 -c ./data/P53.cls \
#              -o test/gsea_reprot_2 \
#              -v --no-plot \
#              -t phenotype
Explanation: 4.3 Command line usage
You may also want to use gsea in command line
End of explanation
# txt, gct file input
ss = gp.ssgsea(data="./tests/data/testSet_rand1200.gct",
               gene_sets="./tests/data/randomSets.gmt", 
               outdir='test/ssgsea_report', 
               sample_norm_method='rank', # choose 'custom' for your own rank list
               permutation_num=0, # skip permutation procedure, because you don't need it
               no_plot=True, # skip plotting, because you don't need these figures
               processes=4, format='png', seed=9)
ss.res2d.sort_index().head()
# or assign a dataframe, or Series to ssgsea()
ssdf = pd.read_csv("./tests/data/temp.txt", header=None, sep="\t")
ssdf.head()
# dataframe with one column is also supported by ssGSEA or Prerank
# But you have to set gene_names as index
ssdf2 = ssdf.set_index(0)
ssdf2.head()
type(ssdf2)
ssSeries = ssdf2.squeeze()
type(ssSeries)
# reuse data
df = pd.read_csv("./tests/data/P53_resampling_data.txt", sep="\t")
df.head()
# Series, DataFrame Example
# supports dataframe and series
ssgs = []
for i, dat in enumerate([ssdf, ssdf2, ssSeries, df]):
    sstemp = gp.ssgsea(data=dat, 
                       gene_sets="./tests/data/genes.gmt", 
                       outdir='test/ssgsea_report_'+str(i),
                       scale=False, # set scale to False to get real original ES
                       permutation_num=0, # skip permutation procedure, because you don't need it
                       no_plot=True, # skip plotting, because you don't need these figures
                       processes=4, seed=10,
                       format='png')
    ssgs.append(sstemp)
Explanation: 5. Single Sample GSEA example
Note: When you run ssGSEA, all genes names in your gene_sets file should be found in your expression table
What's ssGSEA? Which one should I use? Prerank or ssGSEA
see FAQ here
5.1 Input format
Assign ssgsea() with a txt file, gct file, pd.DataFrame, or pd.Seires(gene name as index)
End of explanation
# normalized es save to res2d attri
# one sample input
# NES
ssgs[0].res2d.sort_index().head()
Explanation: 5.2 Access Enrichment Score (ES) and NES
results save to two attribute:
obj.resultsOnSamples: ES
obj.res2d: NES
End of explanation
# ES
# convert dict to DataFrame
es = pd.DataFrame(ssgs[-1].resultsOnSamples)
es.sort_index().head()
# if set scale to True, then 
# Scaled ES equal to es/gene_numbers
ses = es/df.shape[0]
ses
# NES 
# scale or no have no affects on final nes value
nes = ssgs[-1].res2d
nes.sort_index().head()
Explanation: Note: 
If you want to obtain the real original enrichment score,
you have to set scale=False
End of explanation
# set --no-scale to obtain the real original enrichment score
# !gseapy ssgsea -d ./data/testSet_rand1200.gct \
#                -g data/temp.gmt \
#                -o test/ssgsea_report2  \
#                -p 4 --no-plot --no-scale
Explanation: 3) command line usage of single sample gsea
End of explanation
# run command inside python console
rep = gp.replot(indir="./tests/data", outdir="test/replot_test")
Explanation: 6. Replot Example
6.1 locate your directory
notes: replot module need to find edb folder to work properly.
keep the file tree like this:
```
data
 |--- edb
 |    |--- C1OE.cls
 |    |--- gene_sets.gmt
 |    |--- gsea_data.gsea_data.rnk
 |    |--- results.edb
```
End of explanation
# !gseapy replot -i data -o test/replot_test
Explanation: 6.2 command line usage of replot
End of explanation
<END_TASK>