引入clusterProfiler与注释数据
library(clusterProfiler)
library(DOSE)
library(DO.db)
library(org.Mm.eg.db)
GO(gene ontology)分析
GO，Gene Ontology，是基因功能国际标准分类体系。它旨在建立一个适用于各种物种的，对基因和蛋白质功能进行限定和描述的，并能随着研究不断深入而更新的语言词汇标准。GO分为分子功能（Molecular Function）、生物过程（Biological Process）、和细胞组成（Cellular Component）三个部分。
ID转换
org.Mm.eg.db数据对象里面包含着各大主流数据库的数据，一般人都比较熟悉的entrez ID 和ensembl 数据库的ID。
#查看ID类型
keytypes(org.Mm.eg.db)
gene <- row.names(diff_gene_deseq2)
#使用select函数进行表格拼接，这是来自数据库的函数
#keys是原始的ID，columns是转换之后的ID，keytype是要指定的原始ID类型
tansid <- select(org.Mm.eg.db,keys = gene,columns = c("SYMBOL","ENTREZID","GENENAME"),keytype = "ENSEMBL")
可视化
#GO分析及画图
ego <- enrichGO(gene = tansid$ENSEMBL,
keyType = "ENSEMBL",
OrgDb = "org.Mm.eg.db",
ont = "BP"
)
# 气泡图
dotplot(ego, font.size=10)
#网络图,但是找不到函数
enrichMap(ego, vertex.label.cex=1.2, layout=igraph::layout.kamada.kawai)
#GO图
plotGOgraph(ego)
横轴为GeneRatio, 代表该GO term下的差异基因个数占差异基因总数的比例，纵轴为富集到的GO Terms的描述信息
有向无环。矩形代表富集到的top10个GO terms, 颜色从黄色过滤到红色，对应p值从大到小。椭圆形代表未富集到前十的GOterms。
##GSEA
# 按照log2FC从大到小排序基因列表
genelist <- diff_gene_deseq2$log2FoldChange
names(genelist) <- rownames(diff_gene_deseq2)
genelist <- sort(genelist, decreasing = TRUE)
#gsea
gsemf <- gseGO(genelist,
OrgDb = org.Mm.eg.db,
keyType = "ENSEMBL",
ont="MF"
)
head(gsemf)
# 画出GSEA图
gseaplot(gsemf, geneSetID="GO:0000977")
GSEA的输入是一个基因表达量矩阵，其中的样本分成了A和B两组，首先对所有基因进行排序，在之前的文章中也有提到排序的标准，这里简单理解就是foldchange, 用来表示基因在两组间表达量的变化趋势。排序之后的基因列表其顶部可以看做是上调的差异基因，其底部是下调的差异基因。
GSEA分析的是一个基因集下的所有基因是否在这个排序列表的顶部或者底部富集，如果在顶部富集，我们可以说，从总体上看，该基因集是上调趋势，反之，如果在底部富集，则是下调趋势。
- RANK IN GENE LIST代表该基因在排序号的列表中的位置
- RANK METRIC SCORE对应排序，比如foldchange值
- RUNNIG ES代表累计的Enrichment score, CORE ENRICHMENT代表是否属于核心基因，即对该基因集的Enerchment score做出了主要贡献的基因。出现基因集中的基因就增加ES值，不出现就减少（上调）。
https://blog.csdn.net/weixin_43569478/article/details/83745105
https://cloud.tencent.com/developer/article/1426130
KEGG
# id转换
x=bitr(rownames(diff_gene_deseq2),fromType = "ENSEMBL",toType = "ENTREZID", OrgDb = "org.Mm.eg.db")
# 获取keggID
kegg <- x[,2]
# KEGG分析，小鼠mmu，http://www.genome.jp/kegg/catalog/org_list.html
ekk <- enrichKEGG(kegg, keyType = "kegg",organism = "mmu", pvalueCutoff = 0.05, pAdjustMethod = "BH", qvalueCutoff = 0.1)
head(summary(ekk))
# 气泡图
barplot(ekk)
cnetplot(ekk, showCategory = 5)
browseKEGG(ekk, "mmu04742")
# 将GO/KEGG结果转换成CSV格式输出
write.csv(as.data.frame(ekk),"KEGG-enrich.csv",row.names =F)
write.csv(as.data.frame(ego),"GO-enrich.csv",row.names =F)
- 横轴为该pathway的差异基因个数，纵轴为富集到的pathway的描述信息
- 灰色的点代表基因，黄色的点代表富集到的pathways, pathways节点的大小对应富集到的基因个数。
- 富集到的差异基因会用红色方框表示,但是没有一个基因比对上去。
https://www.jianshu.com/p/380bc3508d72
