#### 安裝包
install.packages("BiocManager")
BiocManager::install(version = "3.13")
BiocManager::install("gprofiler2")
BiocManager::install("clusterProfiler")
BiocManager::install("AnnotationHub")
BiocManager::install("org.Bt.eg.db")

GO分析（上下調(diào)基因分開做问芬，可用于BP,CC,MF分開畫圖）

## 方法2：下載到本地加載,每次使用上傳,（推薦）
library(AnnotationDbi)
setwd("D:/生信/sheep.db")
sheep.db <- AnnotationDbi::loadDb('org.Ovis_aries.eg.sqlite')
columns(sheep.db)

## 牛
library(AnnotationDbi)
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富集分析,上下調(diào)分開做
#對于加載的注釋庫的使用簿晓，以上述為例症革，就直接在 OrgDb 中指定人（org.Hs.eg.db）或綿羊（sheep）
enrich.go <- enrichGO(gene = genes$V1,  #基因列表文件中的基因名稱
                      OrgDb = "org.Bt.eg.db",  #指定物種的基因數(shù)據(jù)庫
                      keyType = 'SYMBOL',  #指定給定的基因名稱類型，例如這里以 SYMBOL 為例
                      ont = 'CC',  #可選 BP仔掸、MF脆贵、CC，也可以指定 ALL 同時計算 3 者
                      pAdjustMethod = 'BH',  #指定 p 值校正方法
                      pvalueCutoff = 1,  #指定 p 值閾值起暮，不顯著的值將不顯示在結(jié)果中
                      qvalueCutoff = 1  #指定 q 值閾值卖氨，不顯著的值將不顯示在結(jié)果中
                     )
# 例如上述指定 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分析（上下調(diào)一起做，強烈推薦纸厉，可以看整體結(jié)果系吭，并可選取其中的部分畫圖）

#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')
# 轉(zhuǎn)基因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 條形圖（上下調(diào)一起蟆盹，不過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() +  # 翻轉(zhuǎn)坐標
  theme(panel.grid.major.y = element_blank(),panel.grid.minor.y = element_blank())+ #刪除縱向網(wǎng)絡(luò)
  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))+ # 兩側(cè)加上標簽文字
  ggtitle("The most enriched GO-BP terms")  # 添加標題

GO氣泡圖（可上下調(diào)一起做孩灯，也不用將BP,CC,MF分開，最好選擇適當數(shù)量）

從這位老哥那里薅過來的逾滥，這里是原文鏈接(https://blog.csdn.net/sweet_yemen/article/details/125457274)

# 富集氣泡圖
 # BiocManager::install("openxlsx")
library(ggplot2)#這個是一個非常厲害的繪圖R包
library(openxlsx)#這個包是用來導入格式為XLSX的數(shù)據(jù)的

goinput <- read.xlsx("go-barplot.xlsx")#載入GO富集分析的結(jié)果文件
# 這里強調(diào)一下富集后的GeneRatio是以n/n這種形式的峰档，需要換成小數(shù)形式，在excel中使用=--("0 "&A1)
# A1代表對應(yīng)的單元格
# 設(shè)置XY軸
x=goinput$GeneRatio#設(shè)置以哪個數(shù)據(jù)為X軸
y=factor(goinput$Description,levels = goinput$Description)#設(shè)置Y軸
# 繪制畫布
p = ggplot(goinput,aes(x,y))#開始以X軸和Y軸繪制畫布
p#這里可以查看一下繪制的程度了寨昙，學代碼時多看看每個新的對象有助于理解代碼

#繪圖讥巡，大概的意思是，又設(shè)置了除XY軸數(shù)據(jù)（Ratio舔哪，GOterm）之后又設(shè)置了三個維度的數(shù)據(jù)欢顷。
#三個數(shù)據(jù)分別是以count為點的大小，以顏色的由嫩綠到深粉紅代表P值的對數(shù)值捉蚤，以點的性狀代表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, 上下調(diào)一起做（適合看上下調(diào)整體的結(jié)果）

#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')
# 轉(zhuǎn)基因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， 上下調(diào)分開做（適合畫圖）

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) # 轉(zhuǎn)換為entrezID缆巧，KEGG只識別這種ID布持，報錯請檢查自己輸入的文件，對自己的代碼有自信
write.table(covID,'ENTREZID-up.txt',sep = '\t',quote = F) # 寫出去保存,也可以不保存陕悬，看個人需求

#每次打開R計算時题暖，它會自動連接kegg官網(wǎng)獲得最近的物種注釋信息，因此數(shù)據(jù)庫一定都是最新的
search_kegg_organism("taurus", by="scientific_name") # 尋找對應(yīng)動物的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值閾值，不顯著的值將不顯示在結(jié)果中
                          qvalueCutoff = 1) # 指定q值閾值拼岳，不顯著的值將不顯示在結(jié)果中
#輸出結(jié)果
write.table(enrich_kegg@result,'kegg-down.txt',quote = F,sep = '\t')

KEGG富集氣泡圖(上下調(diào)分開做推薦)

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")