對于非模式生物或者無參考基因組的項目廷蓉,經(jīng)常需要進行基因的功能注釋,而GO注釋是基因功能注釋的重要部分马昙。有很多軟件能夠獲得GO注釋的信息桃犬,例如interproscan、eggnog-mapper和blas2go等行楞。
在這里利用interproscan軟甲獲得了水稻的GO注釋信息攒暇,一般需要將注釋信息整理成gene2go,go2gene以及wego的格式子房,有時也要提取GO Level2的GOID進行分析形用。下面是GO注釋結(jié)果的整理及分析過程。
首先证杭,是interproscan的注釋結(jié)果文件和GO.db文件田度,如下:
```
library(tidyverse)
library(GO.db)
# 加載interpro注釋結(jié)果
interpro <- read_tsv("all.pep.tsv", col_names = FALSE)
#從GO.db導出所有GO注釋信息備用
godb <- select(GO.db, keys(GO.db), columns(GO.db))
write_tsv(godb, "godb.txt")
```
得到的GO注釋是以“|”分割的字符串(不同的軟件各不相同,但都是GOID連接起來字符串)躯砰,這種格式不便于后續(xù)分析,如下圖:
另外携丁,從GO.db包中提取出所有的GO信息琢歇,如下圖:
現(xiàn)在將GO注釋結(jié)果整理成gene2go的格式,并和GO.db關(guān)聯(lián)起來梦鉴,
```
# 提取GO注釋結(jié)果生成gene2go格式數(shù)據(jù)(注釋用的是蛋白序列李茫,要將相同的基因合并,具體情況具體分析)
gene2go <- interpro %>% dplyr::select(Gene = X1, GOID = X14) %>%
? mutate(Gene = str_replace(Gene, "\\.\\d+$", "")) %>% na.omit() %>%
? separate(GOID, paste0("X", 1:(max(str_count(.$GOID,"\\|"))+1), seq = ""), sep = "\\|")? %>%
? gather(key = "X", value = "GOID", -Gene) %>% dplyr::select(Gene, GOID) %>%
? na.omit() %>% base::unique()
write_tsv(gene2go, "gene2go.txt")
# gene2go和godb合并生成完整的GO注釋信息
go_annot <- gene2go %>% left_join(godb)
write_tsv(go_annot, "go_annot.txt")
```
gene2go結(jié)果如下圖:
go_annot結(jié)果如下圖:
整理成go2gene格式肥橙,并統(tǒng)計每個GOID的gene數(shù)目魄宏,代碼如下:
```
# 生成go2gene格式數(shù)據(jù)
go2gene <- gene2go %>% group_by(GOID) %>%
? summarise(Gene = str_c(Gene, collapse = ",")) %>%
? mutate(Count = str_count(Gene, ",")+1) %>% arrange(desc(Count)) %>%
? left_join(godb)
write_tsv(go2gene, "go2gene.txt")
```
go2gene結(jié)果如下圖:
生成 wego 格式:
```
wego <- gene2go %>% group_by(Gene) %>%
? summarise(GOID = str_c(GOID, collapse = ",")) %>%
? separate(GOID, paste0("X", 1:(max(str_count(.$GOID,","))+1), seq = ""), sep = ",")
write_tsv(wego, "wego.txt", col_names = FALSE, na = "")
```
結(jié)果如下圖
可以將wego.txt上傳到wego網(wǎng)站進行分析,可以得到分類匯總信息和GO二級分類統(tǒng)計圖存筏,結(jié)果如下
獲取GO Level2的注釋信息宠互,
```
# 獲得GO Level2的go2gene格式注釋信息
# 利用clusterprofile中的getGOLevel函數(shù),該函數(shù)并沒有開放給用戶椭坚,需要自己調(diào)用源代碼
source("~/software/scripts/go-utilities.R")
go_level2 <- union(getGOLevel("MF", 2), getGOLevel("CC", 2)) %>% union(getGOLevel("BP", 2))
go2gene_level2 <- go2gene %>% filter(GOID %in% go_level2)
```
結(jié)果如下:
這里leve2的GOID只有8個予跌,和wego的結(jié)果差別還是很大。