参考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
提交了7万条蛋白序列,注释时间仅用20分钟就完成了,相当厉害。
完成之后,会提供5种格式(csv/xlsx/gff/Orthologs/seeds)的数据下载。都下载下来更好,省的落下了信息。 csv与xlsx中的存储信息是一致的
在后续构建的内容中,使用了csv格式的数据。
在构建之前,需对csv数据做一下简单处理:
(1)去掉文件前几行“##”注释行,
(2)去掉表头最前面的“#”-保留表头;
(3)去掉文件最后面的注释行;
(4)将文件内的空数据“-”,替换为NA。这一步在参考1里没有提及的细节,一定要注意,非常重要。
做好以上的准备工作,就可以就开始构建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运行环境中。
##过程信息如下:
* 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)
>
到此,从蛋白序列注释,到建库就算完成了。