#### 安装包
install.packages("BiocManager")
BiocManager::install(version = "3.13")
BiocManager::install("gprofiler2")
BiocManager::install("clusterProfiler")
BiocManager::install("AnnotationHub")
BiocManager::install("org.Bt.eg.db")
GO分析（上下调基因分开做，可用于BP,CC,MF分开画图）
## 方法2：下载到本地加载,每次使用上传,（推荐）
setwd("D:/生信/sheep.db")
sheep.db <- AnnotationDbi::loadDb('org.Ovis_aries.eg.sqlite')
columns(sheep.db)
## 牛
setwd("D:/生信/cattle.db")
bos.db <- AnnotationDbi::loadDb('org.Bt.eg.db.sqlite')
columns(bos.db)
#读取基因列表文件中的基因名称
genes <- read.table('GENEID.txt', header = T, stringsAsFactors = FALSE,sep = '\t')
library("org.Bt.eg.db")
library(clusterProfiler)
#GO富集分析,上下调分开做
#对于加载的注释库的使用，以上述为例，就直接在 OrgDb 中指定人（org.Hs.eg.db）或绵羊（sheep）
enrich.go <- enrichGO(gene = genes$V1, #基因列表文件中的基因名称
OrgDb = "org.Bt.eg.db", #指定物种的基因数据库
keyType = 'SYMBOL', #指定给定的基因名称类型，例如这里以 SYMBOL 为例
ont = 'CC', #可选 BP、MF、CC，也可以指定 ALL 同时计算 3 者
pAdjustMethod = 'BH', #指定 p 值校正方法
pvalueCutoff = 1, #指定 p 值阈值，不显著的值将不显示在结果中
qvalueCutoff = 1 #指定 q 值阈值，不显著的值将不显示在结果中
)
# 例如上述指定 ALL 同时计算 BP、MF、CC，这里将它们作个拆分后输出
BP <- enrich.go[enrich.go$ONTOLOGY=='BP', ]
CC <- enrich.go[enrich.go$ONTOLOGY=='CC', ]
MF <- enrich.go[enrich.go$ONTOLOGY=='MF', ]
write.table(as.data.frame(BP), 'HUgo.BP.txt', sep = '\t', row.names = FALSE, quote = FALSE)
write.table(as.data.frame(CC), 'HUgo.CC.txt', sep = '\t', row.names = FALSE, quote = FALSE)
write.table(as.data.frame(MF), 'HUgo.MF.txt', sep = '\t', row.names = FALSE, quote = FALSE)
GO分析（上下调一起做，强烈推荐，可以看整体结果，并可选取其中的部分画图）
#up及down基因的合并处理
gene_down <- read.table('up_gene.txt', header = T, stringsAsFactors = FALSE,sep = '\t')
gene_up <- read.table('down_gene.txt', header = T, stringsAsFactors = FALSE,sep = '\t')
# 转基因ID
library(clusterProfiler)
deg.up.gene <- bitr(gene_up$x, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = bos.db)$ENTREZID
deg.down.gene <- bitr(gene_down$x, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = bos.db)$ENTREZID
# 将up 和 down基因合并到一个表格
deg.ei.ls <- list(up = deg.up.gene, down = deg.down.gene)
# GO
library("org.Bt.eg.db")
GO.cmp <- compareCluster(deg.ei.ls,
fun = "enrichGO",
qvalueCutoff = 0.1,
pvalueCutoff = 0.01,
pAdjustMethod = 'BH',
ont = 'BP', #可选 BP、MF、CC，也可以指定 ALL 同时计算 3 者，建议分开做输出，默认按照pvalue升序排列的
OrgDb = "org.Bt.eg.db")
write.table(GO.cmp@compareClusterResult,'GO-up-down.txt',sep = '\t',quote = F)
GO 条形图（上下调一起，不过BP,CC,MF等需要分开做）
library(ggplot2)
library(ggsci)
ddf<-read.table('GO-up-down.txt',header = T,sep = '\t')
dat=ddf
dat$Cluster = ifelse(dat$Cluster =='up',1,-1)
dat$pvalue = -log10(dat$pvalue)
dat$pvalue <- dat$pvalue*dat$Cluster
dat=dat[order(dat$pvalue,decreasing = F),]
ggplot(dat, aes(x=reorder(Description,order(pvalue, decreasing = F)),
y=pvalue,fill=Cluster)) +
geom_bar(stat="identity",width = 0.78,position = position_dodge(0.7)) +
scale_fill_gradient(low="#3685af",high="#af4236",guide = FALSE)+
scale_x_discrete(name ="GO Terms") +
#scale_y_discrete(labels =Description )
scale_y_continuous(name ="-log10Pvalue",expand = c(0,0)) +
coord_flip() + # 翻转坐标
theme(panel.grid.major.y = element_blank(),panel.grid.minor.y = element_blank())+ #删除纵向网络
theme(axis.text.y = element_blank())+ #刻度标签空白
theme(axis.ticks = element_blank())+ #刻度为空白
theme(panel.border = element_blank())+ #边框为空白
theme(axis.text=element_text(face = "bold",size = 15),
axis.title = element_text(face = 'bold',size = 15), #坐标轴字体
plot.title = element_text(size = 20,hjust = 0.3),
legend.position = "top",panel.grid = element_line(colour = 'white'))+
geom_text(aes(label=Description),size=3.5,hjust=ifelse(dat$pvalue>0,1.5,-0.5))+ # 两侧加上标签文字
ggtitle("The most enriched GO-BP terms") # 添加标题
GO气泡图（可上下调一起做，也不用将BP,CC,MF分开，最好选择适当数量）
从这位老哥那里薅过来的，这里是原文链接(https://blog.csdn.net/sweet_yemen/article/details/125457274)
# 富集气泡图
# BiocManager::install("openxlsx")
library(ggplot2)#这个是一个非常厉害的绘图R包
library(openxlsx)#这个包是用来导入格式为XLSX的数据的
goinput <- read.xlsx("go-barplot.xlsx")#载入GO富集分析的结果文件
# 这里强调一下富集后的GeneRatio是以n/n这种形式的，需要换成小数形式，在excel中使用=--("0 "&A1)
# A1代表对应的单元格
# 设置XY轴
x=goinput$GeneRatio#设置以哪个数据为X轴
y=factor(goinput$Description,levels = goinput$Description)#设置Y轴
# 绘制画布
p = ggplot(goinput,aes(x,y))#开始以X轴和Y轴绘制画布
p#这里可以查看一下绘制的程度了，学代码时多看看每个新的对象有助于理解代码
#绘图，大概的意思是，又设置了除XY轴数据（Ratio，GOterm）之后又设置了三个维度的数据。
#三个数据分别是以count为点的大小，以颜色的由嫩绿到深粉红代表P值的对数值，以点的性状代表GO类型（MF,CC,BP）
p1 = p + geom_point(aes(size=Count,color=-0.5*log(pvalue),shape=ONTOLOGY,))+
scale_color_gradient(low = "SpringGreen", high = "DeepPink")
p1
#添加各类标签
p2 = p1 + labs(color=expression(-log[10](pvalue)),
size="Count",
x="Ratio",
y="Go_term",
title="Go enrichment of test Genes")
p2
KEGG, 上下调一起做（适合看上下调整体的结果）
#up及down基因的合并处理
gene_down <- read.table('up_gene.txt', header = T, stringsAsFactors = FALSE,sep = '\t')
gene_up <- read.table('down_gene.txt', header = T, stringsAsFactors = FALSE,sep = '\t')
# 转基因ID
library(clusterProfiler)
deg.up.gene <- bitr(gene_up$x, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = bos.db)$ENTREZID
deg.down.gene <- bitr(gene_down$x, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = bos.db)$ENTREZID
# 将up 和 down基因合并到一个表格
deg.ei.ls <- list(up = deg.up.gene, down = deg.down.gene)
library(clusterProfiler)
R.utils::setOption("clusterProfiler.download.method",'auto') ##一定要用这个，不然报错别怪我哦
kegg.cmp <- compareCluster(deg.ei.ls,
fun = "enrichKEGG",
qvalueCutoff = 1,
pvalueCutoff = 1,
pAdjustMethod = 'BH',
organism = "bta")
dotplot(kegg.cmp)
write.table(kegg.cmp@compareClusterResult,'kegg-up-down.txt',sep = '\t',quote = F)
kegg， 上下调分开做（适合画图）
library(clusterProfiler)
genes <- read.table('txt', header = T, stringsAsFactors = FALSE,sep = '\t')
covID <- bitr(genes$V1, fromType = "SYMBOL",
toType="ENTREZID",OrgDb= bos.db , drop = TRUE) # 转换为entrezID，KEGG只识别这种ID，报错请检查自己输入的文件，对自己的代码有自信
write.table(covID,'ENTREZID-up.txt',sep = '\t',quote = F) # 写出去保存,也可以不保存，看个人需求
#每次打开R计算时，它会自动连接kegg官网获得最近的物种注释信息，因此数据库一定都是最新的
search_kegg_organism("taurus", by="scientific_name") # 寻找对应动物的organism
R.utils::setOption("clusterProfiler.download.method",'auto') ##一定要用这个，不然报错别怪我
enrich_kegg <- enrichKEGG(gene = covID$ENTREZID, # 基因列表文件中的基因名
organism = 'bta', # 指定物种，牛为bta
keyType = 'kegg', # kegg富集
pAdjustMethod = 'BH', # 指定矫正方法
pvalueCutoff = 1, # 指定p值阈值，不显著的值将不显示在结果中
qvalueCutoff = 1) # 指定q值阈值，不显著的值将不显示在结果中
#输出结果
write.table(enrich_kegg@result,'kegg-down.txt',quote = F,sep = '\t')
KEGG富集气泡图(上下调分开做推荐)
library(ggplot2)
dotplot(enrich_kegg,
showCategory = 15,#展示前15个点，默认为10个
color = 'pvalue',
title = "kegg_dotplot")
KEGG富集条形图
barplot(enrich_kegg,
showCategory = 15,#展示前15个点，默认为10个
color = 'pvalue',
title = "kegg_dotplot")
网友评论