R|clusterProfiler-富集分析
简单总结clusterProfiler包进行GO、KEGG的富集分析方法,结果输出及内置的图形展示。
一 bioconductor包安装
#PS:安装是否顺利,看缘分了
#source("https://bioconductor.org/biocLite.R")
#biocLite("DOSE")
#biocLite("topGO")
#biocLite("clusterProfiler")
#biocLite("pathview")
##加载需要的R包
library(DOSE)
library(org.Hs.eg.db)
library(topGO)
library(clusterProfiler)
library(pathview)
二、数据载入及转化
需要将输入的基因格式转为clusterProfiler可分析的格式,通过功能函数bitr实现各种ID之间的转化。通过keytypes函数可查看所有支持及可转化类型
keytypes(org.Hs.eg.db)
1.向量方式读入#适合少量基因
x <- c("AASDH","ABCB11","ADAM12","ADAMTS16","ADAMTS18")
test = bitr(x, #数据集
fromType="SYMBOL", #输入为SYMBOL格式
toType="ENTREZID", # 转为ENTERZID格式
OrgDb="org.Hs.eg.db") #人类 数据库
head(test,2) SYMBOL ENTREZID
1 AASDH 132949
2 ABCB11 8647
2.基因list文件读入
data <- read.table("gene",header=FALSE) #单列基因名文件
data$V1 <- as.character(data$V1) #需要character格式,然后进行ID转化
#将SYMBOL格式转为ENSEMBL和ENTERZID格式
test1 = bitr(data$V1, fromType="SYMBOL", toType=c("ENSEMBL", "ENTREZID"), OrgDb="org.Hs.eg.db")
head(test1,2) SYMBOL ENSEMBL ENTREZID
1 AASDH ENSG00000157426 132949
2 ABCB11 ENSG00000073734 8647
3. 内置示例数据
data(geneList, package="DOSE") #富集分析的背景基因集
gene <- names(geneList)[abs(geneList) > 2]
gene.df <- bitr(gene, fromType = "ENTREZID", toType = c("ENSEMBL", "SYMBOL"), OrgDb = org.Hs.eg.db)
head(gene.df,2) ENTREZID ENSEMBL SYMBOL
1 4312 ENSG00000196611 MMP1
2 8318 ENSG00000093009 CDC45
三、 GO分析
3.1 groupGO 富集分析
ggo <- groupGO(gene = test1$ENTREZID, OrgDb = org.Hs.eg.db, ont = "CC",level = 3,readable = TRUE)
3.2 enrichGO 富集分析
ego_ALL <- enrichGO(gene = test1$ENTREZID,
universe = names(geneList), #背景基因集
OrgDb = org.Hs.eg.db, #没有organism="human",改为OrgDb=org.Hs.eg.db
#keytype = 'ENSEMBL',
ont = "ALL", #也可以是 CC BP MF中的一种
pAdjustMethod = "BH", #矫正方式 holm”, “hochberg”, “hommel”, “bonferroni”, “BH”, “BY”, “fdr”, “none”中的一种
pvalueCutoff = 1, #P值会过滤掉很多,可以全部输出
qvalueCutoff = 1,
readable = TRUE) #Gene ID 转成gene Symbol ,易读
head(ego_ALL,2)
ONTOLOGY ID Description
GO:0002887 BP GO:0002887 negative regulation of myeloid leukocyte mediated immunity
GO:0033004 BP GO:0033004 negative regulation of mast cell activation
GeneRatio BgRatio pvalue p.adjust qvalue geneID Count
GO:0002887 2/121 10/11461 0.004706555 0.7796682 0.7796682 CD300A/CD84 2
GO:0033004 2/121 10/11461 0.004706555 0.7796682 0.7796682 CD300A/CD84 2
其中:ONTOLOGY:CC BP MF
GO ID: Gene Ontology数据库中唯一的标号信息
Description :Gene Ontology功能的描述信息
GeneRatio:差异基因中与该Term相关的基因数与整个差异基因总数的比值
BgRation:所有( bg)基因中与该Term相关的基因数与所有( bg)基因的比值
pvalue: 富集分析统计学显著水平,一般情况下, P-value < 0.05 该功能为富集项
p.adjust 矫正后的P-Value
qvalue:对p值进行统计学检验的q值
geneID:与该Term相关的基因
Count:与该Term相关的基因数
3.2.1 setReadable函数进行转化
ego_MF <- enrichGO(gene = test1$ENTREZID, universe = names(geneList),OrgDb = org.Hs.eg.db,ont = "MF", pAdjustMethod = "BH",pvalueCutoff = 1,qvalueCutoff = 1,readable = FALSE)
ego_MF1 <- setReadable(ego_MF, OrgDb = org.Hs.eg.db)
3.3 GSEA分析 (暂略)
3.4 输出结果及图像
3.4.1 输出结果
write.csv(summary(ego_ALL),"ALL-enrich.csv",row.names =FALSE)
3.4.2 绘制图形
##可视化--点图
dotplot(ego_MF,title="EnrichmentGO_MF_dot")#点图,按富集的数从大到小的
##可视化--条形图
barplot(ego_MF, showCategory=20,title="EnrichmentGO_MF")#条状图,按p从小到大排,绘制前20个Term
##可视化--
plotGOgraph(ego_MF)
3.5 过滤
3.5.1 利用内置cutoff设置阈值,by设定指标
#可能会有问题,还不太清楚
ego_MF.fil <- simplify(ego_MF, cutoff=0.05, by="pvalue", select_fun=min)
3.5.2 存储结果然后过滤
ego_ALL.sig <- ego_ALL[ego_ALL$pvalue <= 0.01]
过滤后为数据框,不能用自带的参数直接绘制,可以使用ggplot2进行绘制。(暂略)
四、KEGG分析
4.1 候选基因进行通路分析
kk <- enrichKEGG(gene = test1$ENTREZID,
organism = 'hsa', #KEGG可以用organism = 'hsa'
pvalueCutoff = 1)
head(kk,2)
ID Description GeneRatio BgRatio
hsa04750 hsa04750 Inflammatory mediator regulation of TRP channels 5/53 97/7387
hsa04020 hsa04020 Calcium signaling pathway 6/53 182/7387
pvalue p.adjust qvalue geneID Count
hsa04750 0.0006135305 0.08589427 0.0807277 40/3556/3708/5608/79054 5
hsa04020 0.0018078040 0.12654628 0.1189345 493/1129/2066/3707/3708/4842 6
4.2 富集结果及图形展示
4.2.1 结果输出文件
write.csv(summary(kk),"KEGG-enrich.csv",row.names =FALSE)
4.2.1 KEGG 气泡图
dotplot(kk,title="Enrichment KEGG_dot")
4.2.2 查看特定通路图
hsa04750 <- pathview(gene.data = geneList,
pathway.id = "hsa04750", #上述结果中的hsa04750通路
species = "hsa",
limit = list(gene=max(abs(geneList)), cpd=1))
五、注释文件、注释库
如果clusterProfiler包没有所需要物种的内置数据库,可以通过自定义注释文件或者自建注释库的方法进行富集分析。
5.1 自定义注释文件
5.1.1 待富集的基因list
data(geneList, package="DOSE")
deg <- names(geneList)[abs(geneList)>2]
5.1.2 读入注释文件:
## downloaded from http://www.disgenet.org/ds/DisGeNET/results/all_gene_disease_associations.tar.gz
gda <- read.delim("all_gene_disease_associations.tsv")
head(gda,1)geneId geneSymbol diseaseId diseaseName score NofPmids NofSnps source
1 1 A1BG C0001418 Adenocarcinoma 0.002732912 1 0 LHGDN
自定义的两个注释文件(重要)
disease2gene=gda[, c("diseaseId", "geneId")]
disease2name=gda[, c("diseaseId", "diseaseName")]
5.1.3 enricher函数进行富集分析
x = enricher(deg, TERM2GENE=disease2gene, TERM2NAME=disease2name)
head(summary(x),1)
ID Description GeneRatio BgRatio pvalue p.adjust
C3642345 C3642345 Luminal A Breast Carcinoma 11/198 110/17074 6.057328e-08 0.0001488891
qvalue geneID Count
C3642345 0.0001287342 9833/55355/3620/891/9122/2167/3169/5304/2625/9/5241 11
5.2 数据库
5.2.1 查看当前支持的物种
Bioconductor - BiocViews
5.2.2 通过AnnotationHub自建OrgDb注释库并富集
library(AnnotationHub)
hub <- AnnotationHub() #较慢
query(hub, "Cricetulus")
Cgriseus <- hub[["AH48061"]]
sample_gene <- sample(keys(Cgriseus), 100)
str(sample_gene)
sample_test <- enrichGO(sample_gene, OrgDb=Cgriseus, pvalueCutoff=1, qvalueCutoff=1)
head(summary(sample_test))
六、参考资料
Statistical analysis and visualization of functional profiles for genes and gene clusters
以上,clusterProfiler包进行GO、KEGG的富集分析方法,结果输出及内置的图形展示。待补充富集结果ggplot2绘制,GSEA等相关分析。
- Java基础-day09-代码题-对象;类;封装
- MySQL replace into的使用细则(r10笔记第48天)
- Win10下用Anaconda安装TensorFlow
- 【Go 语言社区】跨域问题解决方案:jsonP客户端和服务器代码
- 图;代码轻松理解,代理
- 巧用闪回数据库来查看历史数据 (r10笔记第47天)
- 【Go 语言社区】Golang内存分配
- 小白也能懂的手写体识别
- 【Go 语言社区】浅析javascript的间隔调用和延时调用
- 说说JSON和JSONP,也许你会豁然开朗-转
- C++动态链接库
- Java基础-day11-接口;多态案例练习
- 记一次sql server 性能调优,查询从20秒至2秒
- axis2开发webservice(1)
- JavaScript 教程
- JavaScript 编辑工具
- JavaScript 与HTML
- JavaScript 与Java
- JavaScript 数据结构
- JavaScript 基本数据类型
- JavaScript 特殊数据类型
- JavaScript 运算符
- JavaScript typeof 运算符
- JavaScript 表达式
- JavaScript 类型转换
- JavaScript 基本语法
- JavaScript 注释
- Javascript 基本处理流程
- Javascript 选择结构
- Javascript if 语句
- Javascript if 语句的嵌套
- Javascript switch 语句
- Javascript 循环结构
- Javascript 循环结构实例
- Javascript 跳转语句
- Javascript 控制语句总结
- Javascript 函数介绍
- Javascript 函数的定义
- Javascript 函数调用
- Javascript 几种特殊的函数
- JavaScript 内置函数简介
- Javascript eval() 函数
- Javascript isFinite() 函数
- Javascript isNaN() 函数
- parseInt() 与 parseFloat()
- escape() 与 unescape()
- Javascript 字符串介绍
- Javascript length属性
- javascript 字符串函数
- Javascript 日期对象简介
- Javascript 日期对象用途
- Date 对象属性和方法
- Javascript 数组是什么
- Javascript 创建数组
- Javascript 数组赋值与取值
- Javascript 数组属性和方法
- 如何在 Linux 中查找一个命令或进程的执行时间
- Ubuntu 18.04 LTS中配置IP地址的完整步骤
- Linux系统下Nginx支持ipv6配置的方法
- 微信研发体系下的分布式配置系统设计概要
- Linux双网卡绑定脚本的方法示例
- Serverless 有一百种玩法,比好玩更好玩
- 如何在容器服务中获取客户端真实源IP
- Linux服务器间文件实时同步的实现
- centos7 设置grub密码及单用户登录实例代码
- Linux命令行快速技巧之定位一个文件的方法
- Linux低电量自动关机的实现方法
- ubuntn备份方法总结(四种)
- linux让程序开机自动运行最简单的方法
- centos克隆linux虚拟机的完整步骤分享
- CentOS7.4下MySQL5.7.28二进制方式安装的方法步骤