[GSEAPY]用 Python 进行基因组富集分析
前言
在生物信息学数据分析中,许多分析软件都是基于R开发的。这里介绍一个可以在Python 中进行基因富集分析的Python 软件 GSEAPY
(Gene Set Enrichment Analysis in Python)
GSEApy is a python wrapper for GESA and Enrichr.
It’s used for convenient GO enrichments and produce publication-quality figures from python.
GSEAPY
安装
可以通过conda
或 pip
进行安装
# if you have conda
$ conda install -c conda-forge -c bioconda gseapy
# or use pip to install the latest release
$ pip install gseapy
pip 安装要是遇到这样的报错
data = self.read(amt=amt, decode_content=decode_content)
File "/opt/conda/lib/python3.9/site-packages/pip/_vendor/urllib3/response.py", line 541, in read
raise IncompleteRead(self._fp_bytes_read, self.length_remaining)
File "/opt/conda/lib/python3.9/contextlib.py", line 135, in __exit__
self.gen.throw(type, value, traceback)
File "/opt/conda/lib/python3.9/site-packages/pip/_vendor/urllib3/response.py", line 443, in _error_catcher
raise ReadTimeoutError(self._pool, None, "Read timed out.")
pip._vendor.urllib3.exceptions.ReadTimeoutError: HTTPSConnectionPool(host='files.pythonhosted.org', port=443): Read timed out.
可以使用清华镜像,进行安装:
$ pip install gseapy -i https://pypi.tuna.tsinghua.edu.cn/simple
富集分析
背景信息
- gene set, 指一组具有相同特征的基因。如一个GO term 对应的多个基因,一个kegg pathway对应的多个基因
- gene set library,多个相关的gene set 。如所有GO term组成一个gene set library.
- enrichment analysis, gene set library 作为注释基因集合,已知的先验知识。对于一个输入基因集合,富集分析通过计算分析哪些注释gene set 显著存在于输入基因集合中。例如:GO 富集分析中,查看哪些GO terms 显著存在于输入基因列表中。
有多种基因集富集分析策略,我们常说的GO/KEGG 富集分析 应该大多数指over represent analysis(ORA)。还有个常用富集方法叫GSEA(Gene Set Enrichment Analysis), 翻译过来也是基因集富集分析。下文GSEA,特指这种策略。
ORA
测试数据,可以从GSEApy/tests/data下载。
富集的函数是enricher
.
先展示一下,富集的代码:
gene_list="./gene_list.txt"
gene_sets='KEGG_2016'
gene_sets=['KEGG_2016','KEGG_2013']
enr = gp.enrichr(gene_list=gene_list,
gene_sets=gene_sets,
organism='Human', # don't forget to set organism to the one you desired! e.g. Yeast
description='kegg',
outdir='test/enrichr',
# no_plot=True,
cutoff=0.5 # test dataset, use lower value from range(0,1)
)
运行完后,'test/enrichr'
目录下存放着会有富集的图片以及文本。
(base) jovyan@95c3096ad9ae:~$ ll test/enrichr
-rw-r--r-- 1 jovyan users 22003 Dec 26 14:59 KEGG_2013.Human.enrichr.reports.pdf
-rw-r--r-- 1 jovyan users 22130 Dec 26 14:59 KEGG_2013.Human.enrichr.reports.txt
-rw-r--r-- 1 jovyan users 25722 Dec 26 14:59 KEGG_2016.Human.enrichr.reports.pdf
-rw-r--r-- 1 jovyan users 48458 Dec 26 14:59 KEGG_2016.Human.enrichr.reports.txt
查看KEGG_2016.Human.enrichr.reports.pdf,图片只显示了前10个,这是由参数top_term=10,所决定的
同时富集也结果保存在enr.results
里,如查看前五个数据
enr.results.head(5)
输出
Gene_set Term Overlap P-value Adjusted P-value Old P-value Old Adjusted P-value Odds Ratio Combined Score Genes
0 KEGG_2016 Osteoclast differentiation Homo sapiens hsa04380 28/132 3.104504e-13 7.885440e-11 0 0 6.659625 191.802220 LILRA6;ITGB3;LILRA2;LILRA5;PPP3R1;FCGR3B;SIRPA...
1 KEGG_2016 Tuberculosis Homo sapiens hsa05152 31/178 4.288559e-12 5.446470e-10 0 0 5.224941 136.763196 RAB5B;ITGB2;PPP3R1;HLA-DMA;FCGR3B;HLA-DMB;CASP...
2 KEGG_2016 Phagosome Homo sapiens hsa04145 28/154 1.614009e-11 1.366528e-09 0 0 5.490501 136.437381 ATP6V1A;RAB5B;ITGB5;ITGB3;ITGB2;HLA-DMA;FCGR3B...
3 KEGG_2016 Rheumatoid arthritis Homo sapiens hsa05323 19/90 2.197884e-09 1.395656e-07 0 0 6.554453 130.668081 ATP6V1A;ATP6V1G1;ATP6V0B;TGFB1;ITGB2;FOS;ITGAL...
4 KEGG_2016 Leishmaniasis Homo sapiens hsa05140 17/73 3.132614e-09 1.591368e-07 0 0 7.422186 145.336773 TGFB1;IFNGR1;PRKCB;IFNGR2;ITGB2;FOS;MAPK14;HLA...
查看enricher函数帮助文档
help(gp.enrichr)
Help on function enrichr in module gseapy.enrichr:
enrichr(gene_list, gene_sets, organism='human',
description='', outdir='Enrichr',
background='hsapiens_gene_ensembl',
cutoff=0.05, format='pdf', figsize=(8, 6),
top_term=10, no_plot=False, verbose=False)
......
......
由帮助文档可知enricher
函数所需参数如下:
- gene_list, 所需查询gene_list,可以是一个列表,也可为文件(一列,每行一个基因)
- gene_sets, gene set library。该参数,有两种形式:
- 可以设置enricher自带的gene set library 详细列表可见https://maayanlab.cloud/Enrichr/#libraries。可单个'KEGG_2016',或多个['KEGG_2016','KEGG_2013']
- 一种自定义gene set library。可以是gmt文件,或者输入一个字典
gene_sets={'term_A':['gene1', 'gene2',...],
'term_B':['gene2', 'gene4',...], ...}
- organism,支持(human, mouse, yeast, fly, fish, worm), 自定义gene_set 则无影响。
- description,工作运行描述
- outdir; 输出目录
- background: 背景基因
- 可以是一个背景基因列表
- 或者一个背景基因数目
- 又或者Biomart dataset name.
- cutoff; pvalue阈值
- format, 输出图片格式('pdf','png','eps'...)
- figsize, 图片大小, (width,height). Default: (6.5,6).
- no_plot:是否不做图
绘图
gseapy 也提供了绘图函数进行绘制
# 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',)
enr.res2d
存储着最近一次查询富集的结果, 上面的例子中, enr.res2d
储存的是'KEGG_2013']富集结果,因为它是list最后一个.
gene_sets=['KEGG_2016','KEGG_2013']
enr.results
有着所有的富集结果,所以我么也可以挑选数据可视化
barplot(enr.results.loc[enr.results["Gene_set"] == "KEGG_2016",], title='KEGG_2016',)
气泡图也是有的;
GSEA
Prerank
Prerank
用于已经排好序的数据来做GSEA。如,根据logFC 从大到小排好序后,去做GSEA。
# gsea_data.gsea_data.rnk 是已经排好序的数据
rnk = pd.read_csv("./gsea_data.gsea_data.rnk", header=None, sep="\t")
rnk.head()
0 | 1 |
---|---|
CTLA2B | 2.502482 |
SCARA3 | 2.095578 |
LOC100044683 | 1.116398 |
pre_res = gp.prerank(rnk=rnk, gene_sets='KEGG_2016',
processes=4,
outdir='test/prerank_report_kegg', format='png', seed=6)
pre_res.res2d.head()
绘图
from gseapy.plot import gseaplot
terms = pre_res.res2d.index
# 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]])
未完待续...
参考
https://gseapy.readthedocs.io/en/latest/introduction.html
推荐阅读
-
【视频】Python用LSTM长短期记忆神经网络对不稳定降雨量时间序列进行预测分析|数据分享|附代码数据
-
Python多元回归实战:用SPSS进行多元回归分析
-
实战教程:用Python从雅虎财经获取数据并进行分析与可视化
-
行动中的数据分析--用 Python 对博客评论数据进行情感分析
-
用 Python 进行词云可视化,带您分析《海贼王》、《火影忍者》和《死神》这三部经典动漫
-
包婷婷 (201550484)作业一 统计软件简介与数据操作-SPSS(Statistical Product and Service Solutions),"统计产品与服务解决方案"软件。最初软件全称为"(SolutionsStatistical Package for the Social Sciences),但是随着SPSS产品服务领域的扩大和服务深度的增加,SPSS公司已于2000年正式将英文全称更改为"统计产品与服务解决方案",标志着SPSS的战略方向正在做出重大调整。为IBM公司推出的一系列用于统计学分析运算、数据挖掘、预测分析和决策支持任务的软件产品及相关服务的总称SPSS,有Windows和Mac OS X等版本。 1984年SPSS总部首先推出了世界上第一个统计分析软件微机版本SPSS/PC+,开创了SPSS微机系列产品的开发方向,极大地扩充了它的应用范围,并使其能很快地应用于自然科学、技术科学、社会科学的各个领域。世界上许多有影响的报刊杂志纷纷就SPSS的自动统计绘图、数据的深入分析、使用方便、功能齐全等方面给予了高度的评价。 R统计软件介绍 R是一套完整的数据处理、计算和制图软件系统。其功能包括:数据存储和处理系统;数组运算工具(其向量、矩阵运算方面功能尤其强大);完整连贯的统计分析工具;优秀的统计制图功能;简便而强大的编程语言:可操纵数据的输入和输出,可实现分支、循环,用户可自定义功能。 与其说R是一种统计软件,还不如说R是一种数学计算的环境,因为R并不是仅仅提供若干统计程序、使用者只需指定数据库和若干参数便可进行一个统计分析。R的思想是:它可以提供一些集成的统计工具,但更大量的是它提供各种数学计算、统计计算的函数,从而使使用者能灵活机动的进行数据分析,甚至创造出符合需要的新的统计计算方法。 该语言的语法表面上类似 C,但在语义上是函数设计语言(functional programming language)的变种并且和Lisp 以及 APL有很强的兼容性。特别的是,它允许在"语言上计算"(computing on the language)。这使得它可以把表达式作为函数的输入参数,而这种做法对统计模拟和绘图非常有用。 R是一个免费的*软件,它有UNIX、LINUX、MacOS和WINDOWS版本,都是可以免费下载和使用的。在R主页那儿可以下载到R的安装程序、各种外挂程序和文档。在R的安装程序中只包含了8个基础模块,其他外在模块可以通过CRAN获得。 二、R语言 R是用于统计分析、绘图的语言和操作环境。R是属于GNU系统的一个*、免费、源代码开放的软件,它是一个用于统计计算和统计制图的优秀工具。 R作为一种统计分析软件,是集统计分析与图形显示于一体的。它可以运行于UNIX,Windows和Macintosh的操作系统上,而且嵌入了一个非常方便实用的帮助系统,相比于其他统计分析软件,R还有以下特点: 1.R是*软件。这意味着它是完全免费,开放源代码的。可以在它的网站及其镜像中下载任何有关的安装程序、源代码、程序包及其源代码、文档资料。标准的安装文件身自身就带有许多模块和内嵌统计函数,安装好后可以直接实现许多常用的统计功能。[2] 2.R是一种可编程的语言。作为一个开放的统计编程环境,语法通俗易懂,很容易学会和掌握语言的语法。而且学会之后,我们可以编制自己的函数来扩展现有的语言。这也就是为什么它的更新速度比一般统计软件,如,SPSS,SAS等快得多。大多数最新的统计方法和技术都可以在R中直接得到。[2] 3. 所有R的函数和数据集是保存在程序包里面的。只有当一个包被载入时,它的内容才可以被访问。一些常用、基本的程序包已经被收入了标准安装文件中,随着新的统计分析方法的出现,标准安装文件中所包含的程序包也随着版本的更新而不断变化。在另外版安装文件中,已经包含的程序包有:base一R的基础模块、mle一极大似然估计模块、ts一时间序列分析模块、mva一多元统计分析模块、survival一生存分析模块等等.[2] 4.R具有很强的互动性。除了图形输出是在另外的窗口处,它的输入输出窗口都是在同一个窗口进行的,输入语法中如果出现错误会马上在窗口口中得到提示,对以前输入过的命令有记忆功能,可以随时再现、编辑修改以满足用户的需要。输出的图形可以直接保存为JPG,BMP,PNG等图片格式,还可以直接保存为PDF文件。另外,和其他编程语言和数据库之间有很好的接口。[2] 5.如果加入R的帮助邮件列表一,每天都可能会收到几十份关于R的邮件资讯。可以和全球一流的统计计算方面的专家讨论各种问题,可以说是全世界最大、最前沿的统计学家思维的聚集地.[2] R是基于S语言的一个GNU项目,所以也可以当作S语言的一种实现,通常用S语言编写的代码都可以不作修改的在R环境下运行。 R的语法是来自Scheme。R的使用与S-PLUS有很多类似之处,这两种语言有一定的兼容性。S-PLUS的使用手册,只要稍加修改就可作为R的使用手册。所以有人说:R,是S-PLUS的一个“克隆”。 但是请不要忘了:R是免费的(R is free)。R语言源代码托管在github,具体地址可以看参考资料。[3] 。 R语言的下载可以通过CRAN的镜像来查找。 R语言有域名为.cn的下载地址,有六个,其中两个由Datagurn,由 中国科学技术大学提供的。R语言Windows版,其中由两个下载地点是Datagurn和 USTC提供的。 三、stata Stata 是一套提供其使用者数据分析、数据管理以及绘制专业图表的完整及整合性统计软件。它提供许许多多功能,包含线性混合模型、均衡重复反复及多项式普罗比模式。用Stata绘制的统计图形相当精美。 新版本的STATA采用最具亲和力的窗口接口,使用者自行建立程序时,软件能提供具有直接命令式的语法。Stata提供完整的使用手册,包含统计样本建立、解释、模型与语法、文献等超过一万余页的出版品。 除此之外,Stata软件可以透过网络实时更新每天的最新功能,更可以得知世界各地的使用者对于STATA公司提出的问题与解决之道。使用者也可以透过Stata. Journal获得许许多多的相关讯息以及书籍介绍等。另外一个获取庞大资源的管道就是Statalist,它是一个独立的listserver,每月交替提供使用者超过1000个讯息以及50个程序。 四、PYTHON
-
用 python 进行可靠性分析 python 可以进行可靠性和有效性数据分析
-
高考即将来临!用 Python 抓取历年高考数据并进行分析
-
用 Python 进行数据分析》第 8 章:绘图和可视化笔记
-
[GSEAPY]用 Python 进行基因组富集分析