利用GO/KEGG注释给这些基因赋以“功能标签”
功能注释:查询感兴趣的基因/基因集合参与哪些可能的生命过程,起到了什么作用
1.差异分析筛选基因:MAOA(按照FC排序取top10)(NCBI-GeneID :4128)
2.进入KEGG搜索界面:https://www.genome.jp/kegg/mapper/color.html
3.选择Organism-specific为:hsa
不知道类型的可以选择Taxonomy查询
4.选择Optional use of outside类型为:NCBI-GeneID
5.输入MAOA基因(如格式:4128 red)
可在genecards查询
功能富集分析的原因:
• 一组基因直接注释的结果是得到大量的功能结点。
• 这些功能具有概念上的交叠现象,不利于进一步的精细分析,所以研究人员希望对得到的功能结点加以过滤和筛选,以便获得更有意义的功能信息。
• 富集分析方法通常是分析一组基因在某个功能结点上是否过出现(over-presentation)。由单个基因的注释分析发展到大基因集合的成组分析。和随机
比较,关注的基因集显著注释的功能节点
由于分析的结论是基于一组相关的基因,而不是根据单个基因,所以富集分析方法增加了研究的可靠性,同时也能够识别出与生物现象最相关的生物过程。
rm(list = ls())
options(stringsAsFactors = F)
library(clusterProfiler)
library(org.Hs.eg.db)
library(GSEABase)
library(ggplot2)
library(tidyverse)
packageVersion("clusterProfiler")
## Error in download.KEGG.Path(species)
# https://github.com/YuLab-SMU/clusterProfiler/pull/471
getOption("clusterProfiler.download.method")
#R.utils::setOption("clusterProfiler.download.method",'auto')
options(clusterProfiler.download.method = "wininet")
#options(clusterProfiler.download.method = "wget")
getOption("clusterProfiler.download.method")
# 读取差异分析结果
load(file = "data/Step03-edgeR_nrDEG.Rdata")
ls()
# 提取所有差异表达的基因名
DEG <- DEG_edgeR_symbol[DEG_edgeR_symbol$regulated!="normal", 2]
DEG <- as.character(na.omit(DEG))
head(DEG)
## ===GO数据库, 输出所有结果,后续可根据pvalue挑选结果
ego_CC <- enrichGO(gene=DEG, OrgDb= 'org.Hs.eg.db', keyType='SYMBOL', ont="CC", pvalueCutoff= 1,qvalueCutoff= 1)
ego_MF <- enrichGO(gene=DEG, OrgDb= 'org.Hs.eg.db', keyType='SYMBOL', ont="MF", pvalueCutoff= 1,qvalueCutoff= 1)
ego_BP <- enrichGO(gene=DEG, OrgDb= 'org.Hs.eg.db', keyType='SYMBOL', ont="BP", pvalueCutoff= 1,qvalueCutoff= 1)
p_BP <- barplot(ego_BP,showCategory = 10, label_format = 100) + ggtitle("Biological process")
p_CC <- barplot(ego_CC,showCategory = 10, label_format = 100) + ggtitle("Cellular component")
p_MF <- barplot(ego_MF,showCategory = 10, label_format = 100) + ggtitle("Molecular function")
plotc <- p_BP/p_CC/p_MF
plotc
ggsave('result/6.enrichGO.png', plotc, width = 10,height = 16)
## 使用ggplot2绘制
# 对富集结果使用qvalue从小到大排列,取top10
data <- ego_CC@result %>% top_n(n = 10, wt = -(qvalue))
colnames(data)
p <- ggplot(data = data, aes(x=-log10(qvalue), y=reorder(Description,-log10(qvalue)) )) +
geom_bar(stat = 'identity', fill="OliveDrab4", position="dodge", width =0.5) +
scale_x_continuous(expand = c(0,0)) + # 调整柱子底部与y轴紧贴
xlab(label = "") + ylab(label = "") + # 设置x,y坐标轴标题
theme_classic() + # 设置主题,只保留x,y轴
ggtitle(label = "-LOG(q-value)") + # 设置标题
theme(plot.title = element_text(size = 12,color="black",hjust = 0.5), # 标题居中
axis.title = element_text(size = 15,color ="black"), # 标题字体大小
axis.text = element_text(size= 12,color = "black")) # 坐标轴字体大小
p
ego_BP <- data.frame(ego_BP)
ego_CC <- data.frame(ego_CC)
ego_MF <- data.frame(ego_MF)
write.csv(ego_BP,'result/6.enrichGO_BP.csv')
write.csv(ego_CC,'result/6.enrichGO_CC.csv')
write.csv(ego_MF,'result/6.enrichGO_MF.csv')
## === KEGG
genelist <- bitr(gene=DEG, fromType="SYMBOL", toType="ENTREZID", OrgDb='org.Hs.eg.db')
genelist <- pull(genelist,ENTREZID)
ekegg <- enrichKEGG(gene = genelist, organism = 'hsa', pvalueCutoff = 1, qvalueCutoff = 1)
p1 <- barplot(ekegg, showCategory=10, label_format = 100)
p2 <- dotplot(ekegg, showCategory=10, label_format = 100)
plotc = p1/p2
plotc
ggsave('result/6.enrichKEGG.png', plot = plotc, width = 8, height = 10)
ekegg <- data.frame(ekegg)
write.csv(ekegg,'result/6.enrichKEGG.csv')
## === 其他数据库通路
geneset <- read.gmt("data/MsigDB/v7.4/h.all.v7.4.symbols.gmt")
as.data.frame(table(geneset$term))
geneset$term <- gsub(pattern = "HALLMARK_","", geneset$term)
geneset$term <- str_to_title(geneset$term)
my_path <- enricher(gene=DEG, pvalueCutoff = 1, qvalueCutoff = 1, TERM2GENE=geneset)
p1 <- barplot(my_path, showCategory=15,color = "pvalue", label_format = 100)
p1
ggsave("result/6.enrich_HALLMARK.png", plot = p1, width = 8, height = 7)
my_path <- data.frame(my_path)
write.csv(my_path,"result/6.enrich_HALLMARK.csv")
如果没有筛选到差异表达基因怎么办?
GSEA原理介绍要求输入:
1.全部基因表达谱
2.预先定义的基因集合
# 清空当前环境变量
rm(list = ls())
options(stringsAsFactors = F)
# 加载包
library(GSEABase)
library(clusterProfiler)
library(enrichplot)
library(ggplot2)
library(stats)
# 加载数据
load("data/Step03-edgeR_nrDEG.Rdata")
DEG <- DEG_edgeR_symbol
## 构造GSEA分析数据
# 去掉没有配对上symbol的行
DEG <- DEG[!is.na(DEG$SYMBOL),]
# 去掉重复行
DEG <- DEG[!duplicated(DEG$SYMBOL),]
geneList <- DEG$logFC
names(geneList) <- DEG$SYMBOL
head(geneList)
geneList <- sort(geneList,decreasing = T)
head(geneList)
tail(geneList)
# 选择gmt文件
geneset <- read.gmt("data/MsigDB/h.all.v7.5.1.symbols.gmt")
# 运行,输出全部结果
egmt <- GSEA(geneList, TERM2GENE=geneset, pvalueCutoff = 1, verbose = T)
#出点图
dotplot(egmt,label_format = 100)
#按p值出点图
dotplot(egmt,color="pvalue",label_format = 100)
# 单个通路图
# 按照通路名
gseaplot2(egmt, "HALLMARK_ADIPOGENESIS",
title = "HALLMARK_ADIPOGENESIS")
# 按照行数
gseaplot2(egmt, 5, color="red", pvalue_table = T)
#按第1到第6个出图,不显示p值
gseaplot2(egmt, 1:6, color="red")
# 保存结果
go_gsea <- as.data.frame(egmt@result)
write.csv(go_gsea,"result/6.gsea_go_fc.csv",row.names = F)
原创声明:本文系作者授权腾讯云开发者社区发表,未经许可,不得转载。
如有侵权,请联系 cloudcommunity@tencent.com 删除。
原创声明:本文系作者授权腾讯云开发者社区发表,未经许可,不得转载。
如有侵权,请联系 cloudcommunity@tencent.com 删除。