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

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

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

上面紅框中就是基因的名字，這里，我們已經(jīng)注釋到的基因，形成一個(gè)txt文件，內(nèi)容如下：

1. R包依賴(lài)

下面先載入需要的R包，如果沒(méi)有安裝，需要安裝一下：

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

2. 下載數(shù)據(jù)庫(kù)

到Bioconductor中（https://www./），檢索該物種的數(shù)據(jù)庫(kù)：

常見(jiàn)的物種數(shù)據(jù)庫(kù)如下：直接在Bioconductor中安裝OrgDB的名稱(chēng)就行了。

這里，我們用的是水稻的數(shù)據(jù)庫(kù)，名稱(chēng)為：org.Osativa.eg.db

3. 載入數(shù)據(jù)庫(kù)和讀取基因名文件

「載入數(shù)據(jù)庫(kù)」

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

「讀取基因型文件」

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

4. 將ID匹配GID

將geneID，替換為數(shù)據(jù)庫(kù)中的GID

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

5. 對(duì)基因列表進(jìn)行GO注釋

GO注釋包括：

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

「MF注釋?zhuān)骸?/strong>

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)

結(jié)果文件：

同樣的，CC和BP的GO注釋?zhuān)瑢nt后面的改為CC和BP即可。

「CC的GO注釋?zhuān)骸?/strong>

## 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注釋?zhuān)骸?/strong>

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

其它類(lèi)型的圖：

## 其它類(lèi)型的圖：
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,  #基因列表文件中的基因名稱(chēng)
  keyType = 'kegg',  
  organism = 'dosa', 
  pAdjustMethod = 'fdr',  #指定 p 值校正方法
  pvalueCutoff = 1,  
  qvalueCutoff = 1)  

運(yùn)行日志：

作圖：

barplot(kegg)
dotplot(kegg)