對于非模式生物或者無參考基因組的項目，經(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導(dǎo)出所有GO注釋信息備用

godb <- select(GO.db, keys(GO.db), columns(GO.db))

write_tsv(godb, "godb.txt")

```

interproscan注釋結(jié)果（之和pfam數(shù)據(jù)庫比對）

得到的GO注釋是以“|”分割的字符串（不同的軟件各不相同，但都是GOID連接起來字符串），這種格式不便于后續(xù)分析，如下圖：

GO注釋原始格式

另外，從GO.db包中提取出所有的GO信息，如下圖：

GO.db(展示部分)

現(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é)果如下圖：

gene2go

go_annot結(jié)果如下圖：

GO注釋完整信息（gene2go格式）

整理成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é)果如下圖：

go2gene格式

生成 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格式

可以將wego.txt上傳到wego網(wǎng)站進行分析，可以得到分類匯總信息和GO二級分類統(tǒng)計圖，結(jié)果如下

Summary

GO Level2 bar chart

獲取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é)果如下：

GO Level2的go2gene格式注釋信息

這里leve2的GOID只有8個，和wego的結(jié)果差別還是很大。