大家好，我是鄧飛，上一次，我們介紹如何根據顯著性snp，使用bedtools根據上下游距離，根據gff文件注釋基因。(使用bedtools進行gwas基因注釋)

這一次，介紹一下如何根據注釋的基因，進行富集分析，主要是看一下GWAS定位的基因有沒有某一個趨勢，也算是一種驗證的方法。比如籽粒大小找到的30個候選基因，如果都與籽粒發(fā)育相關的生化途徑一致，那就說明找到的都是相關的基因。

之前用于注釋基因需要的gff文件：

上面紅框中就是基因的名字，這里，我們已經注釋到的基因，形成一個txt文件，內容如下：

1. R包依賴

下面先載入需要的R包，如果沒有安裝，需要安裝一下：

  library(clusterProfiler)
  library(enrichplot)
  library(topGO)
  library(Rgraphviz)
  library(openxlsx)
  library(ggplot2)

2. 下載數據庫

到Bioconductor中（https://www./），檢索該物種的數據庫：

常見的物種數據庫如下：直接在Bioconductor中安裝OrgDB的名稱就行了。

這里，我們用的是水稻的數據庫，名稱為：org.Osativa.eg.db

3. 載入數據庫和讀取基因名文件

「載入數據庫」

library(org.Osativa.eg.db)
db <- org.Osativa.eg.db 
organism <- "dosa" # 物種的名稱

「讀取基因型文件」

geneid = read.csv("gene_total.txt",header = F)
head(geneid)

4. 將ID匹配GID

將geneID，替換為數據庫中的GID

map_id = AnnotationDbi::select(db, keys = geneid, columns=c("GID"), keytype = "RAP")
head(map_id)

5. 對基因列表進行GO注釋

GO注釋包括：

• MF注釋
• CC注釋
• BP注釋

「MF注釋：」

go_MF =enrichGO(map_id$GID, 
                 OrgDb=db,
                 keyType = "GID",
                 ont="MF", 
                 pvalueCutoff=1,
                 qvalueCutoff=1, 
                 pAdjustMethod="none")
write.xlsx(go_MF,"go_MF.xlsx")
dotplot(go_MF,color="pvalue")
ggsave("go_MF_dotplot.pdf",width=12,height=6)

結果文件：

同樣的，CC和BP的GO注釋，將ont后面的改為CC和BP即可。

「CC的GO注釋：」

## CC
go_CC =enrichGO(map_id$GID, 
                OrgDb=db,
                keyType = "GID",
                ont="CC", 
                pvalueCutoff=1,
                qvalueCutoff=1, 
                pAdjustMethod="none")
write.xlsx(go_CC,"go_CC.xlsx")
dotplot(go_CC,color="pvalue")
ggsave("go_CC_dotplot.pdf",width=12,height=6)

「BP的GO注釋：」

## BP
go_BP =enrichGO(map_id$GID, 
                 OrgDb=db,
                 keyType = "GID",
                 ont="BP", 
                 pvalueCutoff=1,
                 qvalueCutoff=1, 
                 pAdjustMethod="none")
write.xlsx(go_BP,"go_BP.xlsx")
dotplot(go_BP,color="pvalue")
ggsave("go_BP_dotplot.pdf",width=12,height=6)

其它類型的圖：

## 其它類型的圖：
barplot(go_BP)
heatplot(go_BP)

6. KEGG富集分析

把基因型的ID后面加上“-01”，并且把g變?yōu)閠

rap_id <- paste0(geneid, "-01")
rap_id <- gsub("g","t",rap_id)
head(rap_id)

富集分析：

geneid = read.csv("a1.txt",header = F)$V1
rap_id <- paste0(geneid, "-01")
rap_id <- gsub("g","t",rap_id)
head(rap_id)
kegg <- enrichKEGG(
  gene = rap_id,  #基因列表文件中的基因名稱
  keyType = 'kegg',  
  organism = 'dosa', 
  pAdjustMethod = 'fdr',  #指定 p 值校正方法
  pvalueCutoff = 1,  
  qvalueCutoff = 1)  

運行日志：

作圖：

barplot(kegg)
dotplot(kegg)