分析 | GO 富集分析

2022-07-16 本文已影响0人 shwzhao

从研一来到组里，一直听到 GO 富集分析几个字。直到现在，研二基本结束了，我都没做过，也不会做。
有一个大概的认识，就是，自己的基因集中某种功能基因的占比要高于这种功能的基因在所有基因中的占比。有很多人解释得更清楚，比如 GO分析学习笔记、徐洲更的基因表达分析（中）- 富集分析、转录组入门(8): 富集分析。

那到底怎么做呢？我看了一些的博客，大体总结了一下：

不用写代码
1.1 在线做，一些公司提供了这样的服务，但可能需要付费；
1.2 TBtools 等软件。

需要写代码，用 R 包 clusterProfiler，但也要分情况
2.1 该物种具有 OrgDb，确定好版本后可以直接下载使用；
2.2 没有 OrgDb，自己构建；
2.3 没有 OrgDb，不用构建。

其他（严谨一点）

我并不研究模式植物，而且已有的 OrgDb 可能存在版本问题。所以 2.2 或 2.3 是待选方法。但是！但是！在是否需要构建以及如何构建 OrgDb 上，我又有了疑问。

在是否需要构建的问题上，我看到徐洲更在功能注释后如何做富集分析中提到“你不需要构建Orgdb，因为Orgdb的用途是进行基因编号和GO/KEGG的转换。你可以直接导入基因号和GO/KEGG编号的对应关系到R里面，然后用clusterProfiler进行数据分析”。

在如何构建的问题上，网上也有许多文章进行了介绍。构建 OrgDb 时，需要gene_info 和 gene2go。gene_info需要有两列信息GID和Gene_name。gene2go需要三列信息GID、GO和EVIDENCE。
在众多博客中，都是用 EggNOG 注释所研究物种的蛋白质序列，再从注释结果中提取信息。我在重复的过程中，发现了其中的不一致，这个Gene_name到底是什么？详细回顾非模式物种注释构建过程中选择的是seed_ortholog列；构建自己物种的orgDb中选择的是eggNOG annot列；使用AnnotationForge包轻松构建非模式物种Orgdb包中选择的是不知道哪里冒出来的X.4列；还有选择Prefered_name ...... 有些列甚至不同版本的 EggNOG 结果可能都不存在/不一致。

我的问题是，Gene_name选择什么重要吗？这个东西是为 OrgDb 的其他功能服务的？如果只是进行 GO 富集分析，这个并不影响？由于我的不求甚解，没有对比结果，也没有深究 OrgDb 到底还能干什么。也许研四快结束的时候，就弄清楚了。还有，因为没有其他证据，大家都把EVIDENCE定义为“IEA”。

那么最终，我选择 2.3 方法做一遍，顺便记录一下很多人使用的 2.2 的方法。
主要参考：使用clusterProfiler对非模式植物进行注释，这篇是在Y叔的公众号上发布的，还提供了一些脚本，非常nice；何同学在勘误-一文掌握转录因子的富集分析（更详细）中只用了TF，应该是同样的道理，都用的enricher()，不同的是一个基因只对应一个TF注释，而可能对应多个GO号。

1. EggNOG

用 EggNOG 进行蛋白功能注释，可以在线注释，也可以下载软件和数据库后本地注释，所需文件只是蛋白质文件。
来了，更好的解释，应该是最好的eggnog-mapper功能注释教程、利用pep序列构建OrgDb进行富集分析（1）。
但是呢，现在版本不一样了，一些操作也不一样了。

1.1 软件及数据下载

软件下载（eggnog-mapper）

git clone https://github.com/eggnogdb/eggnog-mapper
cd eggnog-mapper

命令行下载的参数，可能是网络原因，我下载不了

$ ./download_eggnog_data.py -h
usage: download_eggnog_data.py [-h] [-D] [-P] [-M] [-H] [-d HMMER_DBS] [-y] [-f] [-s] [-q] [--data_dir]

optional arguments:
  -h, --help    show this help message and exit
  -D            Do not install the diamond database (default: False)
  -P            Install the Pfam database, required for de novo annotation or realignment (default: False)
  -M            Install the MMseqs2 database, required for "emapper.py -m mmseqs" (default: False)
  -H            Install the HMMER database specified with "-d TAXID". Required for "emapper.py -m hmmer -d TAXID" (default: False)
  -d HMMER_DBS  Tax ID of eggNOG HMM database to download. e.g. "-H -d 2" for Bacteria. Required if "-H". Available tax IDs can be found at http://eggnog5.embl.de/#/app/downloads. (default: None)
  -y            assume "yes" to all questions (default: False)
  -f            forces download even if the files exist (default: False)
  -s            simulate and print commands. Nothing is downloaded (default: False)
  -q            quiet_mode (default: False)
  --data_dir    Directory to use for DATA_PATH. (default: None)

数据本地下载（http://eggnog5.embl.de/download/）后上传、解压

$ ls ./data
eggnog.db.gz  eggnog_proteins.dmnd.gz  eggnog.taxa.db  eggnog.taxa.db.traverse.pkl  mmseqs.tar.gz  pfam.tar.gz
$ gunzip ./data/*.gz

添加环境变量

$ echo 'PATH=$PATH:'`pwd` >> ~/.bashrc
$ source ~/.bashrc

1.2 蛋白注释

$ /home/shwzhao/miniconda3/bin/python3 emapper.py -h

具体的帮助信息：https://github.com/eggnogdb/eggnog-mapper/wiki/

执行
--cpu: 默认 1
-i: 输入文件，默认蛋白 fasta 序列
-m: 选择注释方式，diamond、mmseqs、hmmer、no_search、cache
--output/-o: 输出文件前缀
Diamond
--dmnd_db: Diamond 数据库文件

python3 emapper.py -m diamond -i pep.fa -o pep.diamond

MMseqs2
--mmseqs_db mmseqs 数据库文件

python3 emapper.py -m mmseqs -i pep.fa -o pep.mmseqs

HMMER
-d: HMMER 数据库前缀，euk、bact、arch

python3 emapper.py -m hmmer -d euk -i pep.fa -o pep.hmmer

2. 实操

2.1 我选择 diamond 方法注释

$ mkdir enrich_analysis
$ cd enrich_analysis
$ python3 ~/bin/software/eggnog-mapper/emapper.py -i ../Arabidopsis_thaliana.pep.fa --cpu 20 -m diamond -o pep.diamond --output_dir .
$ ls -lh
total 26M
-rw-rw-r-- 1 shwzhao shwzhao  18M Jul 15 20:01 pep.diamond.emapper.annotations
-rw-rw-r-- 1 shwzhao shwzhao 6.3M Jul 15 19:59 pep.diamond.emapper.hits
-rw-rw-r-- 1 shwzhao shwzhao 1.8M Jul 15 20:00 pep.diamond.emapper.seed_orthologs

来看一下结果大概长什么样子

$ head -6 pep.diamond.emapper.annotations
## Fri Jul 15 19:59:58 2022
## emapper-2.1.7
## /home/shwzhao/bin/software/eggnog-mapper/emapper.py -i /home/train/public_data/genome/Arabidopsis_thaliana/pep.faa --cpu 20 -m diamond -o pep.diamond --output_dir .
##
#query  seed_ortholog   evalue  score   eggNOG_OGs      max_annot_lvl   COG_category    Description     Preferred_name  GOs     EC      KEGG_ko KEGG_Pathway    KEGG_Module     KEGG_Reaction   KEGG_rclassBRITE   KEGG_TC CAZy    BiGG_Reaction   PFAMs
AT1G01010       72658.Bostr.5325s0009.1.p       4.37e-212       597.0   2D32C@1|root,2SPZ5@2759|Eukaryota,3805H@33090|Viridiplantae,3GQ0B@35493|Streptophyta,3HUB1@3699|Brassicales     35493|StreptophytaSNAC domain-containing protein 69-like   -       -       -       -       -       -       -       -       -       -       -       -       NAM

$ tail -4 pep.diamond.emapper.annotations
ATMG01410       3702.ATMG01410.1        1.78e-144       407.0   2CKHQ@1|root,2S5E8@2759|Eukaryota,37WKE@33090|Viridiplantae,3GKNN@35493|Streptophyta    35493|Streptophyta      S       BEST Arabidopsis thaliana protein match is Mitovirus RNA-dependent RNA polymerase (TAIR    -       -       -       -       -       -       -       -       -       -       -       -       Mitovir_RNA_pol
## 25632 queries scanned
## Total time (seconds): 62.11912512779236
## Rate: 412.63 q/s

2.2 GO 富集分析

下载 go-basic.obo，处理生成 go.tb

$ wget -c http://snapshot.geneontology.org/ontology/go-basic.obo
$ du -sh go-basic.obo
30M     go-basic.obo
$ awk -F ": " 'BEGIN{print "GO\tDescription\tlevel"}{if($0~/\[Term\]/)term=NR;if(NR==term+1)id=$2;if(NR==term+2)name[id]=$2;if(NR==term+3)namespace[id]=$2}END{for(i in name)print i"\t"name[i]"\t"namespace[i]}' go-basic.obo > go.tb

文件准备，解析 eggnog 结果，联合 go.tb

library(tidyverse)
library(clusterProfiler, lib = "/home/zcli/R/x86_64-pc-linux-gnu-library/4.0/")

emapper <- read_tsv("pep.diamond.emapper.annotations", comment = "##") %>% rename(query = "#query")
# 读取时去掉`##`开头的行，而`#query`开头的行作为表头

GOinfo <- read_tsv("go.tb")

GOannotation <- emapper %>%
  dplyr::select(GID = query, GO = GOs) %>%
  dplyr::filter(GO != "-") %>%
  separate_rows(GO, sep = ",") %>%
  left_join(GOinfo, by = "GO") %>%
  dplyr::select(-Description)

GOannotation <- split(GOannotation, GOannotation$level)

提取自己的genelist（注意是一个 vector），然后用 enricher() 富集
这里根据参考，把TERM2GENE根据不同GO类型分开运行。应该是可以同时运行的吧？是的吧？啊？

TF <- read_tsv("Arabidopsis_thaliana.TF.txt", col_names = F)

genelist <- TF %>%
  head(1000) %>%
  .$X1

enrich <- enricher(genelist,
                   TERM2GENE = GOannotation[['molecular_function']][c(2,1)],
                   TERM2NAME = GOinfo[1:2])

绘图

我的建议是 as.data.frame() 或 as_tibble() 后用 ggplot 画。
目前还不知道怎么用网络的形式呈现富集结果。

pdf("dotplot.pdf")
barplot(enrich)
dotplot(enrich)
# 注意有时需要调整字体：font.size
dev.off()

其他
富集分析结果呈现还可以有很多操作。比如 simpleEnrichment--对GO富集分析结果进一步处理！ ......

补充：记录 OrgDb 构建的方式

# 文件准备
library(tidyverse)

emapper <- read_tsv("pep.diamond.emapper.annotations", comment = "##") %>% rename(query = "#query")

gene2go <- emapper %>%
  dplyr::select(GID = query, GO = GOs) %>%
  dplyr::filter(GO != "-") %>%
  separate_rows(GO, sep = ",") %>%
  mutate(EVIDENCE = "IEA")

gene_info <- emapper %>%
  select(GID = query, Gene_name = seed_ortholog)

# 构建OrgDb，首先得到`"./org.Atha.eg.db"`
AnnotationForge::makeOrgPackage(gene_info = gene_info,
                 go = gene2go,
                 maintainer = 'shwzhao <shwzhao997@bjfu.edu.cn>',
                 author = 'shwzhao',
                 outputDir = ".",
                 tax_id = 1842,
                 genus = 'Ar',
                 species = 'tha',
                 goTable = "go",
                 version = "1.0")

# 封装包，得到 `"./org.Atha.eg.db_1.0.tar.gz"`
pkgbuild::build('./org.Atha.eg.db', dest_path = ".")

# 安装包，得到`org.Atha.eg.db`
install.packages('org.Atha.eg.db',
                 repos = NULL,
                 type = "source",
                 lib = "./orgdb")

# 调用包
library(org.Atha.eg.db, lib = "./orgdb")

# enrichGO() 进行富集分析
# `gene`: entrez 基因 id 的向量
# `OrgDb`: 
# `universe`: 背景基因
# `keyType`: 
# `ont`: (`"MF"` | `"BP"` | `"CC"`) 分别为分子功能，生物学过程，细胞组分
# `pAdjustMethod`: (`"holm"` | `"hochberg"` | `"hommel"` | `"bonferroni"` | `"BH"` | `"BY"` | `"fdr"` | `"none"`)
# `pvalueCutoff`: p 值
# `qvalueCutoff`: q 值
# `minGSSize`: 
# `maxGSSize`: 
# `readable`: 是否将基因 id 转换成基因名