GO富集分析
2022-12-01 本文已影响0人
Jason数据分析生信教室
转录组分析传送门
NGS手把手教学之零基础RNA-seq转录组分析实践，两套方案(2022年最新）
通路富集分析简介
GO富集详解 clusterProfiler包
KEGG富集详解(待更新)
Reactome富集详解(待更新)
富集分析结果可视化大全(待更新)
目录
- GO归类
- GO过表现分析(ORA)
- GO 基因集合富集分析(GSEA)
- 非模型生物的GO分析
- 可视化GO
1. GO归类
前文提到GO分三种，molecular function (MF), biological process (BP), and cellular component (CC)。
在clusterProfiler包里，groupGO() 可以用来给基因功能分类。本例中用到了DOSE包里提供的数据集geneList 。
BiocManager::install("clusterProfiler")
BiocManager::install("org.Hs.eg.db")
library(org.Hs.eg.db)
library(clusterProfiler)
data(geneList, package="DOSE")
gene <- names(geneList)[abs(geneList) > 2]
# Entrez gene ID
head(gene)
显示gene ID
## [1] "4312" "8318" "10874" "55143" "55388" "991"
GO归类
ggo <- groupGO(gene = gene,
OrgDb = org.Hs.eg.db,
ont = "CC", #也可以是MF,BP
level = 3,
readable = TRUE)
head(ggo)
## ID Description Count GeneRatio geneID
## GO:0000133 GO:0000133 polarisome 0 0/207
## GO:0000408 GO:0000408 EKC/KEOPS complex 0 0/207
## GO:0000417 GO:0000417 HIR complex 0 0/207
## GO:0000444 GO:0000444 MIS12/MIND type complex 0 0/207
## GO:0000808 GO:0000808 origin recognition complex 0 0/207
## GO:0000930 GO:0000930 gamma-tubulin complex 0 0/207
2. GO过表现分析(ORA)
clusterProfiler包的话用enrichGO()就可以
ego <- enrichGO(gene = gene,
universe = names(geneList),
OrgDb = org.Hs.eg.db,
ont = "CC",
pAdjustMethod = "BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05,
readable = TRUE)
head(ego)
## ID Description GeneRatio
## GO:0005819 GO:0005819 spindle 26/201
## GO:0000779 GO:0000779 condensed chromosome, centromeric region 16/201
## GO:0072686 GO:0072686 mitotic spindle 17/201
## GO:0000775 GO:0000775 chromosome, centromeric region 18/201
## GO:0098687 GO:0098687 chromosomal region 23/201
## GO:0000776 GO:0000776 kinetochore 15/201
## BgRatio pvalue p.adjust qvalue
## GO:0005819 306/11853 1.072029e-11 3.151766e-09 2.888837e-09
## GO:0000779 114/11853 7.709944e-11 8.659125e-09 7.936756e-09
## GO:0072686 133/11853 8.835841e-11 8.659125e-09 7.936756e-09
## GO:0000775 158/11853 1.684987e-10 1.179661e-08 1.081250e-08
## GO:0098687 272/11853 2.006225e-10 1.179661e-08 1.081250e-08
## GO:0000776 106/11853 2.733425e-10 1.339378e-08 1.227644e-08
## geneID
## GO:0005819 CDCA8/CDC20/KIF23/CENPE/ASPM/DLGAP5/SKA1/NUSAP1/TPX2/TACC3/NEK2/CDK1/MAD2L1/KIF18A/BIRC5/KIF11/TRAT1/TTK/AURKB/PRC1/KIFC1/KIF18B/KIF20A/AURKA/CCNB1/KIF4A
## GO:0000779 CENPE/NDC80/HJURP/SKA1/NEK2/CENPM/CENPN/ERCC6L/MAD2L1/KIF18A/CDT1/BIRC5/TTK/NCAPG/AURKB/CCNB1
## GO:0072686 KIF23/CENPE/ASPM/SKA1/NUSAP1/TPX2/TACC3/CDK1/MAD2L1/KIF18A/KIF11/TRAT1/AURKB/PRC1/KIFC1/KIF18B/AURKA
## GO:0000775 CDCA8/CENPE/NDC80/TOP2A/HJURP/SKA1/NEK2/CENPM/CENPN/ERCC6L/MAD2L1/KIF18A/CDT1/BIRC5/TTK/NCAPG/AURKB/CCNB1
## GO:0098687 CDCA8/CENPE/NDC80/TOP2A/HJURP/SKA1/NEK2/CENPM/RAD51AP1/CENPN/CDK1/ERCC6L/MAD2L1/KIF18A/CDT1/BIRC5/EZH2/TTK/NCAPG/AURKB/CHEK1/CCNB1/MCM5
## GO:0000776 CENPE/NDC80/HJURP/SKA1/NEK2/CENPM/CENPN/ERCC6L/MAD2L1/KIF18A/CDT1/BIRC5/TTK/AURKB/CCNB1
## Count
## GO:0005819 26
## GO:0000779 16
## GO:0072686 17
## GO:0000775 18
## GO:0098687 23
## GO:0000776 15
只要是数据库OrgDb 支持的gene ID格式都可以直接使用。但是需要在keyType 里指定gene ID的格式。
gene.df <- bitr(gene, fromType = "ENTREZID",
toType = c("ENSEMBL", "SYMBOL"),
OrgDb = org.Hs.eg.db)
ego2 <- enrichGO(gene = gene.df$ENSEMBL,
OrgDb = org.Hs.eg.db,
keyType = 'ENSEMBL',
ont = "CC",
pAdjustMethod = "BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05)
head(ego2, 3)
## ID Description GeneRatio
## GO:0005819 GO:0005819 spindle 30/233
## GO:0072686 GO:0072686 mitotic spindle 21/233
## GO:0000779 GO:0000779 condensed chromosome, centromeric region 17/233
## BgRatio pvalue p.adjust qvalue
## GO:0005819 422/21916 1.954166e-16 5.413039e-14 4.175744e-14
## GO:0072686 179/21916 3.911063e-16 5.416822e-14 4.178662e-14
## GO:0000779 154/21916 7.525516e-13 6.948560e-11 5.360280e-11
## geneID
## GO:0005819 ENSG00000134690/ENSG00000117399/ENSG00000137807/ENSG00000138778/ENSG00000066279/ENSG00000126787/ENSG00000154839/ENSG00000262634/ENSG00000137804/ENSG00000088325/ENSG00000013810/ENSG00000117650/ENSG00000170312/ENSG00000164109/ENSG00000121621/ENSG00000089685/ENSG00000138160/ENSG00000163519/ENSG00000112742/ENSG00000178999/ENSG00000198901/ENSG00000237649/ENSG00000233450/ENSG00000056678/ENSG00000204197/ENSG00000186185/ENSG00000112984/ENSG00000087586/ENSG00000134057/ENSG00000090889
## GO:0072686 ENSG00000137807/ENSG00000138778/ENSG00000066279/ENSG00000154839/ENSG00000262634/ENSG00000137804/ENSG00000088325/ENSG00000013810/ENSG00000170312/ENSG00000164109/ENSG00000121621/ENSG00000138160/ENSG00000163519/ENSG00000178999/ENSG00000198901/ENSG00000237649/ENSG00000233450/ENSG00000056678/ENSG00000204197/ENSG00000186185/ENSG00000087586
## GO:0000779 ENSG00000138778/ENSG00000080986/ENSG00000123485/ENSG00000154839/ENSG00000262634/ENSG00000117650/ENSG00000100162/ENSG00000166451/ENSG00000186871/ENSG00000164109/ENSG00000121621/ENSG00000167513/ENSG00000089685/ENSG00000112742/ENSG00000109805/ENSG00000178999/ENSG00000134057
## Count
## GO:0005819 30
## GO:0072686 21
## GO:0000779 17
3. GO 基因集合富集分析(Gene Set Enrichment Analysis, GSEA)
gseGO() 可以实现GSEA分析。
ego3 <- gseGO(geneList = geneList,
OrgDb = org.Hs.eg.db,
ont = "CC",
minGSSize = 100,
maxGSSize = 500,
pvalueCutoff = 0.05,
verbose = FALSE)
4. 非模型生物的GO分析
如果需要处理非模型生物也就是没有OrgDb 的时候可以自己通过AnnotationHub搜索到自己需要的数据库。如果实在没有的话可以用biomaRt或者 Blast2GO来获取GO注释。又或者用AnnotationForge 直接创建一个自己的数据库。
5. 可视化GO
p值，关系网络，参与程度都会有显示出来。
goplot(ego)
Fig1.png
