口腔医学小站

口腔文献阅读与资源分享

R包clusterProfiler: 转换ID、GO/KEGG富集分析

2021年3月28日 发表者 kimi1006 | 暂无评论

R包clusterProfiler: 转换ID、GO/KEGG富集分析

12020.04.23 16:21:40字数 1,707阅读 3,773

参考这篇
ID转换不用怕，clusterProfiler帮你忙
这篇详细讲了无参考基因组该怎么办，值得一看

简单介绍一下几种常用的ID：
Ensemble id：一般以ENSG开头，后边跟11位数字。如TP53基因：ENSG00000141510
Entrez id：通常为纯数字。如TP53基因：7157
Symbol id：为我们常在文献中报道的基因名称。如TP53基因的symbol id为TP53
Refseq id：NCBI提供的参考序列数据库：可以是NG、NM、NP开头，代表基因，转录本和蛋白质。如TP53基因的某个转录本信息可为NM_000546

1 转换ID

1.1 载入包library(clusterProfiler)

library(clusterProfiler)
library(DOSE)
library(org.Hs.eg.db)
246

org.Hs.eg.db是人类基因组注释，如果需要其他物种的注释信息，可以去bioconductor这里获得

1.2. keytypes()查看注释包中支持的ID转换类型，基本包括了所以常用的ID类型

keytypes()

1.3. bitr()函数转换ID

函数全部内容如下

bitr(geneID, fromType, toType, OrgDb, drop = TRUE
2

geneID：输入待转换的geneID
fromType：输入的ID类型
toType：希望输出的ID类型
OrgDb：注释对象的信息
drop：drop = FALSE 保留空值

以TP53为例，希望输出’ENTREZID’,’ENSEMBL’,’REFSEQ’

my_id <- c("TP53")

output <- bitr(my_id,
               fromType = 'SYMBOL',
               toType = c('ENTREZID','ENSEMBL','REFSEQ'),
               OrgDb = 'org.Hs.eg.db')
24681012
转换后的结果

1.4. bitr_kegg()函数进行基因ID与蛋白质ID的转换

函数全部内容如下，与bitr()类似

bitr_kegg(geneID, fromType, toType, organism, drop = TRUE)
2

在参数这里与bitr()有些区别！
geneID：输入待转换的geneID
fromType：输入的ID类型只能是kegg（与entrez id相同）, ncbi-geneid, ncbi-proteinid or uniprot中的一个
toType：希望输出的ID类型，也只能是kegg（与entrez id相同）, ncbi-geneid, ncbi-proteinid or uniprot中的一个
organism：hsa，代表人类，其他的物种名称可以参考kegg的网站
drop：去除空值与否

TP53的entrez id为7157

kegg <- bitr_kegg(7157,
                   fromType = 'kegg',
                   toType = c('ncbi-geneid', 'ncbi-proteinid', 'uniprot'),
                   organism = 'hsa')
2468

报错了，所以一次应该只能进行一种转换！

改成这样↓

kegg <- bitr_kegg(7157,
                  fromType = 'kegg',
                  toType = 'ncbi-proteinid',
                  organism = 'hsa')

246810

转换为uniprot

kegg <- bitr_kegg(7157,
                  fromType = 'kegg',
                  toType = 'uniprot',
                  organism = 'hsa')
2468

出现了三个uniprot

去UniProt网站上查询一下为什么会有三个：
P04637显示的是TP53基因的蛋白质表达水平，级别是Reviewed，就是其来源为UniProtKB/Swiss-Prot。

P04637

K7PPA8显示的是转录本水平的表达，级别是Unreviewed，就是其来源为UniProtKB/TrEMBL

K7PPA8

Q53GA5显示的是转录本水平的表达，级别是Unreviewed，就是其来源为UniProtKB/TrEMBL

Q53GA5

一般ID转换仅仅为开始的准备工作，将自己的数剧转换好之后可以进行后续的分析。
另外，利用clusterProfiler包可以进行许多丰富的下游分析，比如GO分析、KEGG分析等等

2 富集分析

参考这篇
最常用的基因注释工具是GO和KEGG注释，这基本上是和差异基因分析一样，必做的事，GO可在BP(生物过程），MF(分子功能），CC(细胞组分）三方面进行注释

GO中的注释信息

KEGG数据库全部的物种及其简写（三个字母）

有很多在线（DAVID、Gather、GOrilla,revigo）或者客户端版的工具（IPA（IPA不是用的GO和KEGG数据库）和FUNRICH）可以用来分析差异基因，我主要实践一下clusterProfiler

clusterProfiler提供有两种富集方式：

① ORA（Over-Representation Analysis）差异基因富集分析（不需要表达值，只需要gene name）
② FCS，它的代表就是GSEA，即基因集(gene set)富集分析（不管有无差异，需要全部genes表达值）

① ORA 差异基因富集分析（不需要表达值，只需要gene name）

1) 先读进来差异基因
sig.gene <- read.delim("你的路径/final_result.txt（之前计算出的差异基因文件）",  sep = 't', stringsAsFactors = FALSE, check.names = FALSE)
gene<-sig.gene[,1]
head(gene)

#转换ID
gene.df<-bitr(gene, fromType = "ENSEMBL", 
              toType = c("SYMBOL","ENTREZID"),
              OrgDb = org.Hs.eg.db,
              drop = FALSE)

head(gene.df)
246810121416182022
2) GO富集
ego_cc<-enrichGO(gene = gene.df$ENSEMBL,
                 OrgDb  = org.Hs.eg.db,
                 keyType = 'ENSEMBL',
                 ont  = "CC",
                 pAdjustMethod = "BH",
                 pvalueCutoff = 0.01,
                 qvalueCutoff = 0.05)

ego_bp<-enrichGO(gene = gene.df$ENSEMBL,
                 OrgDb = org.Hs.eg.db,
                 keyType = 'ENSEMBL',
                 ont  = "BP",
                 pAdjustMethod = "BH",
                 pvalueCutoff = 0.01,
                 qvalueCutoff = 0.05)
24681012141618202224262830

参数意义：
keyType 指定基因ID的类型，默认为ENTREZID
OrgDb 指定该物种对应的org包的名字
ont 代表GO的3大类别，BP, CC, MF，也可是全部ALL

pAdjustMethod 指定多重假设检验矫正的方法，有“ holm”, “hochberg”, “hommel”, “bonferroni”, “BH”, “BY”, “fdr”, “none”中的一种
cufoff 指定对应的阈值
readable=TRUE 代表将基因ID转换为gene symbol。

barplot() #富集柱形图
dotplot() #富集气泡图
cnetplot() #网络图展示富集功能和基因的包含关系
emapplot() #网络图展示各富集功能之间共有基因关系
heatplot() #热图展示富集功能和基因的包含关系

3) barplot()作图显示，依赖于ggplot2包
library(ggplot2)
barplot(ego_bp,showCategory = 18,title="The GO_BP enrichment analysis of all DEGs ")  
24

ggplot2参数意义具体可以看这篇

barplot()

dotplot()

emapplot()
4) 或者是KEGG分析

输入的gene要为vector格式的entrezid，可以用is.vector()判断是否是vector格式

library(stringr)

kk<-enrichKEGG(gene = gene.df$ENTREZID,
               organism = "hsa",
               keyType = "kegg",          
               pAdjustMethod = "BH",
               pvalueCutoff = 0.01)
2468101214

KEGG

但是这样产生的kk几乎都是Na，更别提画图了
于是我去看了Y叔的手册，发现需要使用DOSE构建，改成

library(DOSE)

data(geneList, package = "DOSE")
gene <- names(geneList)[abs(geneList) > 2]

kk<-enrichKEGG(gene = gene,
               organism = "hsa",
               keyType = "kegg",
               pAdjustMethod = "BH",
               pvalueCutoff = 0.01)
2468101214161820

就可以了，继续往下进行

#bar图
barplot(kk,showCategory = 25, title="The KEGG enrichment analysis of all DEGs")+
  scale_size(range=c(2, 12))+
  scale_x_discrete(labels=function(kk) str_wrap(kk,width = 25))
#点图
dotplot(kk,showCategory = 25, title="The KEGG enrichment analysis of all DEGs")+
  scale_size(range=c(2, 12))+
  scale_x_discrete(labels=function(kk) str_wrap(kk,width = 25))
246810121416
barplot()

dotplot()

cnetplot()

heatplot()

② GSEA 基因集(gene set)富集分析（不管有无差异，需要全部genes表达值）

这篇为网页版操作
好处：可以发现被差异基因舍弃的genes可能参与了某重要生理过程或信号通路
要点：用所有差异基因来跑GSEA，而不是按筛选条件筛选出来的，否则GSEA就失去了意义

1）读入全部差异基因数据
library(dplyr)
all_deseq <- read.delim("你的路径/DESeq2-res.txt",  sep = 't', stringsAsFactors = FALSE, check.names = FALSE)
24
部分内容
2） 提取出gene_id 和 log2FoldChange列
ensemblid_log2fdc <-dplyr::select(all_deseq, gene_id, log2FoldChange)
2
3) 可以看到这时gene_id还不是整数形式，所以对这列gsub()取整，命名为symbol列添加到ensemblid_log2fdc中
ensemblid_log2fdc$symbol <- gsub("\.\d*", "", all_deseq[,2])
2
4）转换id：将ENSEMBLID转为ENTREZID
entrezid <-bitr(ensemblid_log2fdc$symbol, 
              fromType = "ENSEMBL", 
              toType = c("ENTREZID"),
              OrgDb = org.Hs.eg.db,
              drop = TRUE)
246810
5）因为6）中merge需要相同列名，所以将symbol列列名改为ENSEMBL
names(ensemblid_log2fdc) <- c("gene_id", "log2FoldChange", "ENSEMBL")
2
6）ensemblid_log2fdc 和 entrezid拼接到一起，依靠的是它们具有相同列名的列——ENSEMBL
final <- merge(ensemblid_log2fdc, entrezid)
2
7) 提取entrezid 和log2fdc列
library(dplyr)    
entrezid_log2fdc <- dplyr::select(final, ENTREZID, log2FoldChange)
24
8) 用arrange()按log2fdc排序，desc()表示降序
final.sort <- arrange(entrezid_log2fdc, desc(log2FoldChange))
2
9) 构建一下要GSEA分析的数据
genelist <- final.sort$log2FoldChange
names(genelist) <- final.sort[,1]
24
10) 进行GSEA分析，需要entrezid，因此我前面才会进行ID转化
gsemf <- gseGO(genelist,
           OrgDb = org.Hs.eg.db,
           keyType = "ENTREZID",
           ont="BP",
           pvalueCutoff = 1)
head(gsemf)
24681012
11）画出GSEA图
gseaplot(gsemf, 1)
2
GSEA图

3 可视化

1.GO富集分析结果可视化

barplot

barplot(ego, showCategory = 10) #默认展示显著富集的top10个，即p.adjust最小的10个

dotplot

dotplot(ego, showCategory = 10)

DAG有向无环图

plotGOgraph(ego) #矩形代表富集到的top10个GO terms, 颜色从黄色过滤到红色，对应p值从大到小。

igraph布局的DAG

goplot(ego)

GO terms关系网络图（通过差异基因关联）

emapplot(ego, showCategory = 30)

GO term与差异基因关系网络图

cnetplot(ego, showCategory = 5)

2.Pathway富集分析结果可视化

barplot

barplot(kk, showCategory = 10)

dotplot

dotplot(kk, showCategory = 10)

pathway关系网络图（通过差异基因关联）

emapplot(kk, showCategory = 30)

pathway与差异基因关系网络图

cnetplot(kk, showCategory = 5)

pathway映射

browseKEGG(kk, “hsa04934”) #在pathway通路图上标记富集到的基因，会链接到KEGG官网

禁止转载，如需转载请通过简信或评论联系作者。

25人点赞

更多精彩内容下载简书APP
“小礼物走一走，来简书关注我”

还没有人赞赏，支持一下

总资产2 (约0.15元)共写了1.3W字获得98个赞共87个粉丝

点击数：0

发表回复

*为必填字段！