更新于2021.4.9
又写了一个script，用来一次性构建多个物种包的。
成功倒是成功了，就是可能效率比较低。
主要是用了for循环，然后用列表分别存储各自物种的数据的。
可能还会有更好的方法，以后再琢磨琢磨
R刚入门，最近在学习非模式物种注释包构建，记录一下。主要参考的教程如下：
相关版本信息:
| 包 | 版本号 |
|---|---|
| R | 4.0.3 |
| Bioconductor | 3.12 |
| AnnotationForge | 1.32.0 |
| AnnotationDbi | 1.52.0 |
| ClusterProfiler | 3.18.1 |
| dplyr | 1.0.3 |
| stringr | 1.4.0 |
| data.table | 1.14.0 |
ppps:我在安装clusterProfiler的时候，出现了安装失败的问题。
如果你也碰到了可以试试修改一下镜像哈~
Tools->Global Options->packages
或者直接用命令设置也可以~
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
##clusterProfiler&AnnotationForge
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("clusterProfiler")
BiocManager::install("AnnotationForge")
##other packages you may don't have
if(! require("dplyr")) install.packages("dplyr")
if(! require("stringr")) install.packages("stringr")
if(! require("data.table")) install.packages("data.table")
方法一：AnnotationHub
https://bioconductor.org/packages/release/bioc/html/AnnotationHub.html
不过这个方法我没有去做，直接用下面的方法自己构建的。想了解的话可以康康上面参考的教程。
方法二：AnnotationForge
- 获取蛋白序列，eggnog-mapper人工构建注释。
这一步我直接用的导师给我的数据，所以没做。可以康康上面参考的教程。
为了方便看，我把教程的代码copy过来。
##下载你想要的序列，教程里的例子是芝麻Sesame
wget http://www.sesame-bioinfo.org/SesameFG/BLAST_search/G608_contig_2014-08-29.FgeneSH.pep.rar
# 解压后上传
# 前提是自己安装好eggnog-mapper并且下载好相应的数据库
emapper.py -m diamond \
-i sesame.fa \
-o diamond \
--cpu 19
# 得到如下信息，然后进行处理，只保留表头query_name这一行的注释信息，去掉头尾的# 等信息
sed -i '/^# /d' diamond.emapper.annotations
sed -i 's/#//' diamond.emapper.annotations
- 加载包&读入文件
library(clusterProfiler)
library(dplyr)
library(stringr)
library(data.table)
options(stringsAsFactors = F)
##这个命令hin重要！！千万不要忘记设置！！
##这样设置，所有在数据框中的字符会被当做character来处理。
#设置工作路径
setwd("D:/xxxx")#记得修改
#读取文件
#这个UTEX2973是我挑出来的一个物种，记得修改！
data <- fread('UTEX2973.tab',data.table = T)
data[data == ""]<- #将空行替换成NA，方便后续去除
我的文件跟用eggnog-mapper得到的细节有些不太一样，截张图参考一下~我一共有20列。然后我挑了一个物种出来先试着构建一个物种包。
data
- 挑选Locus_tag&Gene
gene_info <- data %>% dplyr::select(GID = Locus_tag,GENENAME = Gene) %>% na.omit()
gene_info
4.Locus_tag&GO, 并得到对应关系
#Step4-1:挑出两列
gterms <- data %>% dplyr::select(GID = Locus_tag,GO_LIST = GO_LIST) %>% na.omit
#Step4-2::得到两者对应关系
#(一)用for循环，效率较低
##先构建一个空数据框GID=》Locus_tag，GO_LIST=》GO号,EVIDENCE =》默认IDA
gene2go <- data.frame(GID = character(),
GO = character(),
EVIDENCE = character())
##然后向其中填充：以GO号为标准，每一行只能有一个GO号，Locus_tag和EVIDENCE可以重复
for(row in 1:nrow(gterms)){
gene_terms <- str_split(gterms[row,"GO_LIST"]," ",simplify = FALSE)[[1]]
gene_id <- gterms[row,"GID"][[1]]
tmp <- data_frame(GID = rep(gene_id,length(gene_terms)),
GO = gene_terms,
EVIDENCE = rep("IEA",length(gene_terms)))
gene2go <- rbind(gene2go,tmp)
}
#（二）用sapply，效率较高
all_go_list = str_split(gterms$GO_LIST," ")
gene2go <- data.frame(GID = rep(gterms$GID,
times = sapply(all_go_list,length)),
GO = unlist(all_go_list),
EVIDENCE = "IEA")
gterms
gene2go
-
挑出Locus_tag&KEGG,得到pathway2name,ko2pathway
这一步需要下载一个json文件
json下载
#Step5-1，挑出Locus_tag&KEGG
gene2ko <- data %>% dplyr::select(GID = Locus_tag,KO = KEGG) %>% na.omit
#Step5-2，pathway2name,ko2pathway
if(!file.exists('kegg_info.RData')){
# 需要下载这个json文件 (经常更新)
# https://www.genome.jp/kegg-bin/get_htext?ko00001
library(jsonlite)
library(purrr)
library(RCurl)
update_kegg <- function(json = "ko00001.json",file = NULL) {
pathway2name <- data.frame(Pathway = character(), Name = character())
ko2pathway <- data.frame(Ko = character(), Pathway = character())
kegg <- fromJSON(json)
for (a in seq_along(kegg[["children"]][["children"]])) {
A <- kegg[["children"]][["name"]][[a]]
for (b in seq_along(kegg[["children"]][["children"]][[a]][["children"]])) {
B <- kegg[["children"]][["children"]][[a]][["name"]][[b]]
for (c in seq_along(kegg[["children"]][["children"]][[a]][["children"]][[b]][["children"]])) {
pathway_info <- kegg[["children"]][["children"]][[a]][["children"]][[b]][["name"]][[c]]
pathway_id <- str_match(pathway_info, "ko[0-9]{5}")[1]
pathway_name <- str_replace(pathway_info, " \\[PATH:ko[0-9]{5}\\]", "") %>% str_replace("[0-9]{5} ", "")
pathway2name <- rbind(pathway2name, tibble(Pathway = pathway_id, Name = pathway_name))
kos_info <- kegg[["children"]][["children"]][[a]][["children"]][[b]][["children"]][[c]][["name"]]
kos <- str_match(kos_info, "K[0-9]*")[,1]
ko2pathway <- rbind(ko2pathway, tibble(Ko = kos, Pathway = rep(pathway_id, length(kos))))
}
}
}
save(pathway2name, ko2pathway, file = file)
}
update_kegg(json = "ko00001.json",file="kegg_info.RData")
}
这一步出过好多问题，查资料查了好久。上面那个超级复杂的代码是参考刘小泽老师的，中途遇到过kegg_info.RData文件没生成出来的问题，修改了一下，然后就得到文件啦！
-
利用gene2ko与ko2pathway将基因与pathway对应起来
gene2ko
ko2pathway
load(file = "kegg_info.RData")
#运行数据框合并前需要做到两个数据框列名对应，将原来的gene2ko中ko修改一下
colnames(ko2pathway)=c("KO",'Pathway')
gene2ko$KO=str_replace(gene2ko$KO,"ko:"," ")
##因为要对应嘛，所以还需要将ko2pathway中的KO列中的K删掉
ko2pathway$KO <- str_replace(ko2pathway$KO, "K", "")
##上面这个处理不是完全一致的，仅仅只是针对我自己的数据。如果你们参考我的这个来做注释包构建的话，建议先观察自己的数据，视情况而定哈~
#合并
gene2pathway <- gene2ko %>% left_join(ko2pathway,by = "KO") %>%
dplyr::select(GID,Pathway) %>% na.omit()
gene2pathway
现在一个个的对应都建立起来啦，可以开始构建注释包啦~
- 构建注释包,AnnotationForge
library(AnnotationForge)
##tax_id可以去网站搜的哈
#https://www.ncbi.nlm.nih.gov/taxonomy/
tax_id="1350461"
genus="Synechococcus"
species="elongatus"
makeOrgPackage(gene_info=gene_info,
go=gene2go,
ko=gene2ko,
maintainer='<xxxxxxxxx@qq.com>',##这里记住<>这俩括号一定要记得加上哈！下面那个也是的！
author='<xxxxxxxxx@qq.com>',
pathway=gene2pathway,
version="0.0.1",
outputDir = ".",
tax_id=tax_id,
genus=genus,
species=species,
goTable="go")
至此，注释包就构建完成啦！！后面的GO、KEGG富集分析啥的我暂时先不写了，现在要继续去完成导师任务了qaq
我用这个代码跑下来是没出啥问题的，如果有问题可以一起讨论哈~
网友评论