使用AnnotationForge包自定义Orgdb包

参考1：https://www.jianshu.com/p/68f7b663c5ff
参考2：https://www.dzbioinformatics.com/2019/09/21/%e6%9e%84%e5%bb%ba%e8%87%aa%e5%ae%9a%e4%b9%89%e6%b3%a8%e9%87%8a%e5%8c%85make-annotation-package/

写在前面的话：
最近在使用clusterProfiler进行GO或者KEGG富集分析时，找不到自己需要的org..eg.db包。这个包会提供每个基因与GO/KEGG的对应关系。
现在R语言找到的org..eg.db包很少。只能自己想办法自定义构建了。
下面的内容是在参考引用网页的基础上做的内容。
参考1虽然很详细了，但自己真照做的时候，有些细节会出问题。因此这里再梳理一下。为后续人避坑。

测试环境：
windows10-64位系统，
R4.0.5
rtools40-x86_64
RstudioVersion 1.2.1335

第一步：通过eggNOG-mapper-5.0网站注释：可以直接提交蛋白序列、核酸序列等，url: http://eggnog-mapper.embl.de/
这里如果提交蛋白序列或核酸序列一定要注意格式。

>a
TTTTTAAAGGG
>b
TTTTTTTTTTTT

避免多余的信息,如果不信，就测试测试吧。
>a    keyword I'm here
TTTTTAAAGGG
>b   keyword I'm here
TTTTTTTTTTTT

image.png image.png image.png

提交了7万条蛋白序列，注释时间仅用20分钟就完成了，相当厉害。
完成之后，会提供5种格式（csv/xlsx/gff/Orthologs/seeds）的数据下载。都下载下来更好，省的落下了信息。 csv与xlsx中的存储信息是一致的

image.png

在后续构建的内容中，使用了csv格式的数据。

在构建之前，需对csv数据做一下简单处理：
（1）去掉文件前几行“##”注释行，
（2）去掉表头最前面的“#”-保留表头；
（3）去掉文件最后面的注释行；
（4）将文件内的空数据“-”，替换为NA。这一步在参考1里没有提及的细节，一定要注意，非常重要。

image.png image.png

做好以上的准备工作，就可以就开始构建orgDB的具体操作了。
第一步：读取数据

setwd("D:/")
##名称由makeOrgPackage函数规定请不要更改
emapper <- read.delim("D:/MM_xbmdxk2v.emapper.annotations.xls") %>%
  dplyr::select(GID=query,Gene_Symbol=Preferred_name, 
                GO=GOs,KO=KEGG_ko,Pathway =KEGG_Pathway, 
                OG=eggNOG_OGs,Gene_Name =seed_ortholog)
#GID=query,这里的query就是csv文件中的表头信息

第二步：从矩阵中提取各种信息，为建库做准备。

#这里共提取了gene_info,  gene2go,gene2ko,gene2pathway,gene2symbol
#你用哪些信息就可以进行相应的增减。
#如果你只是做go富集分析，其实gene_info和gene2go就足够了
#gene2symbol在这里感觉纯粹是凑数，这是参考引文2弄的。
#gene2ko，gene2pathway是用来做kegg富集分析的。

#提取GID与Gene_Name信息，参考1是将X.4作为Gene_Name，应该是自己定义的信息，这里将seed_ortholog作为Gene_Name信息，你可以根据实际再调整。
gene_info <- dplyr::select(emapper,GID,Gene_Name) %>%
  dplyr::filter(!is.na(Gene_Name))

#提取GID与GO信息，组成goTable，建库时的goTable需要三列（GID, GO和EVIDENCE），少一列，就会出错.
gene2go <- dplyr::select(emapper,GID,GO) %>%
  separate_rows(GO, sep = ',', convert = F) %>%
  dplyr::filter(GO!="",!is.na(GO))%>%   #这是只提取有GO注释信息的行，判断的标准时GO信息不是NA，这也就是为什么前面必须将“-”替换为NA，不替换就无法进行有效过滤。
  mutate(EVIDENCE = 'A')     #硬生生加了1列EVIDENCE，全部赋值A,凑数的。
dim(gene2go)    #查看数据维度。
#[1] 1523399       3

#提取GID与KO信息，这里只有2列信息
gene2ko<- dplyr::select(emapper,GID,KO) %>%
 separate_rows(KO, sep = ',', convert = F) %>%
  dplyr::filter(!is.na(KO))
dim(gene2ko)
#[1] 30530     2

#提取GID与Pathway信息，这里只有2列信息
gene2pathway<- dplyr::select(emapper,GID,Pathway) %>%
separate_rows(Pathway, sep = ',', convert = F) %>%
  dplyr::filter(!is.na(Pathway))
 dim(gene2pathway)
#[1] 143056      2

#提取GID与Gene_Symbol信息，Gene_Symbol是Preferred_name信息，这里只有2列信息
gene2symbol<- dplyr::select(emapper,GID,Gene_Symbol) %>%
  dplyr::filter(!is.na(Gene_Symbol))
dim(gene2symbol)
#[1] 4561    2

第三步：如果提取信息没有任何输出错误信息。就可以建库了。

library(AnnotationForge)  #如果没有就安装一下。
AnnotationForge::makeOrgPackage(gene_info=gene_info,
                                go=gene2go,
                                ko=gene2ko,
                                pathway=gene2pathway,
                                symbol=gene2symbol,
                                maintainer='CJT<chaojiangtao@caas.cn>',
                                author='CJT',
                                version="0.1" ,
                                outputDir=".", 
                                tax_id="4097",
                                genus="N",
                                species="t",
                                goTable="go")
#下面输出的过程信息，如有问题就会报错~~~~~~~~~~~~~~~~~~
Populating genes table:
genes table filled
Populating gene_info table:
gene_info table filled
Populating go table:
go table filled
Populating ko table:
ko table filled
Populating pathway table:
pathway table filled
Populating symbol table:
symbol table filled
table metadata filled
'select()' returned many:1 mapping between keys
and columns
Dropping GO IDs that are too new for the current GO.db
Populating go table:
go table filled
Populating go_bp table:
go_bp table filled
Populating go_cc table:
go_cc table filled
Populating go_mf table:
go_mf table filled
'select()' returned many:1 mapping between keys
and columns
Populating go_bp_all table:
go_bp_all table filled
Populating go_cc_all table:
go_cc_all table filled
Populating go_mf_all table:
go_mf_all table filled
Populating go_all table:
go_all table filled
Creating package in ./org.Nt.eg.db 
Now deleting temporary database file
[1] "./org.Nt.eg.db"
There were 50 or more warnings (use warnings() to see the first 50)
#虽然有警告信息，但还好不影响后续建库。

第四步：对生成的org.MA.eg.db进行封装。
这一步是在DOS命令窗口下运行，运行命令时，转到org.Nt.eg.db所在的目录下运行，或者提供准确而路径也行。

D:\>R CMD build org.Nt.eg.db
###输出的过程信息如下：
* checking for file 'org.Nt.eg.db/DESCRIPTION' ... OK
* preparing 'org.Nt.eg.db':
* checking DESCRIPTION meta-information ... OK
* checking for LF line-endings in source and make files and shell scripts
* checking for empty or unneeded directories
* building 'org.Nt.eg.db_0.1.tar.gz'

到这里就完成建库了。
第五步：将此库导入R运行环境中。

image.png

##过程信息如下：
* installing *source* package 'org.Nt.eg.db' ...
** using staged installation
** R
** inst
** byte-compile and prepare package for lazy loading
** help
*** installing help indices
  converting help for package 'org.Nt.eg.db'
    finding HTML links ... 好了
    org.Nt.egBASE                           html  
    org.Nt.egORGANISM                       html  
    org.Nt.eg_dbconn                        html  
** building package indices
** testing if installed package can be loaded from temporary location
*** arch - i386
*** arch - x64
** testing if installed package can be loaded from final location
*** arch - i386
*** arch - x64
** testing if installed package keeps a record of temporary installation path
* DONE (org.Nt.eg.db)

#尝试导入该库，如果没错误信息，就可以正常使用了。
>library(org.Nt.eg.db)
>

到此，从蛋白序列注释，到建库就算完成了。