常用的分析方法：
- 过表征分析 (over representation analysis, ORA) ：先找出差异基因，再对差异基因进行如GO/KEGG等富集
- 基因富集分析 (gene set enrichment analysis, GSEA)：对所有的基因进行富集，综合考虑每一个基因对于基因集的贡献
GO是基因本体论联合会建立的一个数据库，旨在建立一个适用于各种物种的、对基因和蛋白功能进行限定和描述的，并能随着研究不断深入而更新的语义词汇标准。GO注释分为三大类，分别是：分子生物学功能（Molecular Function，MF）、生物学过程（Biological Process，BP）和细胞学组分（Cellular Components，CC），通过这三个功能大类，对一个基因的功能进行多方面的限定和描述。
KEGG是一个整合了基因组、化学和系统功能信息的综合数据库。KEGG下属4个大类和17和子数据库，其中KEGG Pathway数据库最为常用。
1 GO Analysis
clusterProfiler 包中的 enrichGO() 可以对差异基因进行GO分析
- 流程：差异基因 -- > 基因ID转换 -- > GO富集 -- > 结果可视化
- 需要文件：DGE分析后筛选出的基因
1.1 GO-基因ID的转换
常用ID：
- Ensemble：欧洲生物信息数据库提供 ENSG+11位数字
- Entrez：NCBI提供，纯数字
- Symbol：常用基因名
- Refseq：NCBI提供的参考序列数据，以NG/NM等开头
建议使用 ENTREZID 格式的ID进行GO分析，可通过bitr()函数转换基因ID
bitr(geneID, fromType="", toType="", OrgDb=, drop=TRUE)
geneID: 需要转换的基因ID文件；fromType：geneID文件中的ID类型；toType：输出的ID类型； OrgDb：注释对象； Drop：是否去除空值
e.g. 以DOSE包中的数据geneList 为例
> data(geneList, package="DOSE")
> head(geneList)
4312 8318 10874 55143 55388 991
4.572613 4.514594 4.418218 4.144075 3.876258 3.677857
#选取差异基因
> gene <- names(geneList)[abs(geneList) >2]
# 转换ID
> geneID <- bitr(gene, fromType = "ENTREZID",
+ toType = c("ENSEMBL", "SYMBOL"),
+ OrgDb = org.Hs.eg.db)
> head(geneID)
ENTREZID ENSEMBL SYMBOL
1 4312 ENSG00000196611 MMP1
2 8318 ENSG00000093009 CDC45
3 10874 ENSG00000109255 NMU
4 55143 ENSG00000134690 CDCA8
5 55388 ENSG00000065328 MCM10
6 991 ENSG00000117399 CDC20
# keytypes(org.Hs.eg.db) 可查询人类基因注释包中支持的ID类型
1.2 GO-GO富集分析
enrichGO()函数对差异基因进行GO富集分析
> gene.GO <- enrichGO(gene = geneID$ENTREZID,
OrgDb = org.Hs.eg.db,
keyType = "ENTREZID",
ont = "CC",
pAdjustMethod = "BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05,
readable = T)
> gene.GO[1:2,]
ID Description GeneRatio BgRatio
GO:0005819 GO:0005819 spindle 25/200 347/19717
GO:0072686 GO:0072686 mitotic spindle 15/200 109/19717
pvalue p.adjust qvalue
GO:0005819 1.612727e-14 4.951071e-12 4.498659e-12
GO:0072686 3.213034e-13 4.932007e-11 4.481337e-11
geneID
GO:0005819 CDCA8/CDC20/KIF23/CENPE/ASPM/DLGAP5/SKA1/NUSAP1/TPX2/TACC3/NEK2/CDK1/MAD2L1/KIF18A/BIRC5/KIF11/TTK/AURKB/PRC1/KIFC1/KIF18B/KIF20A/AURKA/CCNB1/KIF4A
GO:0072686 KIF23/CENPE/ASPM/SKA1/NUSAP1/TPX2/TACC3/CDK1/MAD2L1/KIF18A/KIF11/AURKB/KIFC1/KIF18B/AURKA
Count
GO:0005819 25
GO:0072686 15
# 可保存生成的 @result 文件
> write.table(gene.GO@result, file = "geneGO.txt",
sep = "\t", row.names = T,col.names = T)
- 参数：
gene差异基因；keyType指定基因ID的类型，默认为ENTREZID；OrgDb指定该物种对应的org包的名字；ont代表GO的3大类别，BP, CC, MF;pAdjustMethod指定多重假设检验矫正的方法；cufoff指定对应的阈值；readable=TRUE将基因ID转换为gene symbol；universe背景基因，若不指定，可使用数据库中所有基因作为背景
2 KEGG Analysis
流程与GO类似，将函数 enrichGO() 改为 enrichKEGG() 即可
使用参数organism 指定物种，需要查找物种对应的 kegg_code
e.g. 以DOSE包中的数据geneList 为例
> search_kegg_organism("human", by="common_name")
kegg_code scientific_name common_name
1 hsa Homo sapiens human
115 phu Pediculus humanus corporis human body louse
> gene.KEGG <- enrichKEGG(gene = geneID$ENTREZID,
organism = "hsa",
keyType = "kegg",
pAdjustMethod = "BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05 )
3 结果可视化
使用enrichplot 包进行富集结果可视化分析
可视化参考文章
- barplot / dotplot
- heatplot
- cnetplot
- emapplot
- goplot / plotGOgraph
3.1 barplot / dotplot
library(enrichplot)
p1 <- barplot(gene.GO, showCategory = 10)
p2 <- dotplot(gene.GO, showCategory = 10)
plot_grid(p1,p2,ncol = 2)
3.2 heatplot
heatplot(gene.GO, foldChange = geneList)
3.3 cnetplot
显示GO基因集中所富集到的基因，则将该基因与GO集连线
cnetplot(gene.GO, showCategory = 5)
3.4 emapplot
显示富集到的GO基因集之间的关系。若不同基因集见存在重叠基因，则基因集间将会连线
emapplot(gene.GO, showCategory = 20)
3.5 goplot / plotGOgraph
显示GO的有向无环图
goplot(gene.ID)
