GO注释结果整理

对于非模式生物或者无参考基因组的项目，经常需要进行基因的功能注释，而GO注释是基因功能注释的重要部分。有很多软件能够获得GO注释的信息，例如interproscan、eggnog-mapper和blas2go等。

在这里利用interproscan软甲获得了水稻的GO注释信息，一般需要将注释信息整理成gene2go，go2gene以及wego的格式，有时也要提取GO Level2的GOID进行分析。下面是GO注释结果的整理及分析过程。

首先，是interproscan的注释结果文件和GO.db文件，如下：

```

library(tidyverse)

library(GO.db)

# 加载interpro注释结果

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

```

interproscan注释结果（之和pfam数据库比对）

得到的GO注释是以“|”分割的字符串（不同的软件各不相同，但都是GOID连接起来字符串），这种格式不便于后续分析，如下图：

GO注释原始格式

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

GO.db(展示部分)

现在将GO注释结果整理成gene2go的格式，并和GO.db关联起来，

```

# 提取GO注释结果生成gene2go格式数据（注释用的是蛋白序列，要将相同的基因合并，具体情况具体分析）

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结果如下图：

gene2go

go_annot结果如下图：

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

整理成go2gene格式，并统计每个GOID的gene数目，代码如下：

```

# 生成go2gene格式数据

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结果如下图：

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

```

结果如下图

wego格式

可以将wego.txt上传到wego网站进行分析，可以得到分类汇总信息和GO二级分类统计图，结果如下

Summary GO Level2 bar chart

获取GO Level2的注释信息，

```

# 获得GO Level2的go2gene格式注释信息

# 利用clusterprofile中的getGOLevel函数，该函数并没有开放给用户，需要自己调用源代码

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)

```

结果如下：

GO Level2的go2gene格式注释信息

这里leve2的GOID只有8个，和wego的结果差别还是很大。