更新于2021.4.9
又寫了一個script蝌戒，用來一次性構(gòu)建多個物種包的召烂。
成功倒是成功了碱工，就是可能效率比較低。
主要是用了for循環(huán)奏夫，然后用列表分別存儲各自物種的數(shù)據(jù)的怕篷。
可能還會有更好的方法，以后再琢磨琢磨
R剛?cè)腴T酗昼，最近在學(xué)習(xí)非模式物種注釋包構(gòu)建廊谓，記錄一下。主要參考的教程如下：
相關(guān)版本信息:
| 包 | 版本號 |
|---|---|
| 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的時候，出現(xiàn)了安裝失敗的問題呛哟。
如果你也碰到了可以試試修改一下鏡像哈~
Tools->Global Options->packages
或者直接用命令設(shè)置也可以~
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
不過這個方法我沒有去做叠荠，直接用下面的方法自己構(gòu)建的。想了解的話可以康康上面參考的教程竖共。
方法二：AnnotationForge
- 獲取蛋白序列蝙叛，eggnog-mapper人工構(gòu)建注釋。
這一步我直接用的導(dǎo)師給我的數(shù)據(jù)公给，所以沒做借帘。可以康康上面參考的教程淌铐。
為了方便看肺然，我把教程的代碼copy過來。
##下載你想要的序列腿准，教程里的例子是芝麻Sesame
wget http://www.sesame-bioinfo.org/SesameFG/BLAST_search/G608_contig_2014-08-29.FgeneSH.pep.rar
# 解壓后上傳
# 前提是自己安裝好eggnog-mapper并且下載好相應(yīng)的數(shù)據(jù)庫
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重要＝滞校翔！千萬不要忘記設(shè)置！灾前！
##這樣設(shè)置防症，所有在數(shù)據(jù)框中的字符會被當(dāng)做character來處理。
#設(shè)置工作路徑
setwd("D:/xxxx")#記得修改
#讀取文件
#這個UTEX2973是我挑出來的一個物種哎甲，記得修改蔫敲！
data <- fread('UTEX2973.tab',data.table = T)
data[data == ""]<- #將空行替換成NA，方便后續(xù)去除
我的文件跟用eggnog-mapper得到的細節(jié)有些不太一樣炭玫，截張圖參考一下~我一共有20列奈嘿。然后我挑了一個物種出來先試著構(gòu)建一個物種包。
data
- 挑選Locus_tag&Gene
gene_info <- data %>% dplyr::select(GID = Locus_tag,GENENAME = Gene) %>% na.omit()
gene_info
4.Locus_tag&GO, 并得到對應(yīng)關(guān)系
#Step4-1:挑出兩列
gterms <- data %>% dplyr::select(GID = Locus_tag,GO_LIST = GO_LIST) %>% na.omit
#Step4-2::得到兩者對應(yīng)關(guān)系
#(一)用for循環(huán)吞加，效率較低
##先構(gòu)建一個空數(shù)據(jù)框GID=》Locus_tag裙犹，GO_LIST=》GO號,EVIDENCE =》默認IDA
gene2go <- data.frame(GID = character(),
GO = character(),
EVIDENCE = character())
##然后向其中填充：以GO號為標(biāo)準(zhǔn)，每一行只能有一個GO號榴鼎，Locus_tag和EVIDENCE可以重復(fù)
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文件 (經(jīng)常更新)
# 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")
}
這一步出過好多問題，查資料查了好久哩陕。上面那個超級復(fù)雜的代碼是參考劉小澤老師的平项，中途遇到過kegg_info.RData文件沒生成出來的問題，修改了一下悍及，然后就得到文件啦闽瓢！
-
利用gene2ko與ko2pathway將基因與pathway對應(yīng)起來
gene2ko
ko2pathway
load(file = "kegg_info.RData")
#運行數(shù)據(jù)框合并前需要做到兩個數(shù)據(jù)框列名對應(yīng)，將原來的gene2ko中ko修改一下
colnames(ko2pathway)=c("KO",'Pathway')
gene2ko$KO=str_replace(gene2ko$KO,"ko:"," ")
##因為要對應(yīng)嘛心赶，所以還需要將ko2pathway中的KO列中的K刪掉
ko2pathway$KO <- str_replace(ko2pathway$KO, "K", "")
##上面這個處理不是完全一致的扣讼，僅僅只是針對我自己的數(shù)據(jù)。如果你們參考我的這個來做注釋包構(gòu)建的話缨叫，建議先觀察自己的數(shù)據(jù)椭符，視情況而定哈~
#合并
gene2pathway <- gene2ko %>% left_join(ko2pathway,by = "KO") %>%
dplyr::select(GID,Pathway) %>% na.omit()
gene2pathway
現(xiàn)在一個個的對應(yīng)都建立起來啦，可以開始構(gòu)建注釋包啦~
- 構(gòu)建注釋包,AnnotationForge
library(AnnotationForge)
##tax_id可以去網(wǎng)站搜的哈
#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")
至此，注釋包就構(gòu)建完成啦Ｋ龃亍蒸健！后面的GO、KEGG富集分析啥的我暫時先不寫了，現(xiàn)在要繼續(xù)去完成導(dǎo)師任務(wù)了qaq
我用這個代碼跑下來是沒出啥問題的似忧，如果有問題可以一起討論哈~
