NCBI-NT库本地化构建流程

数据库下载

https://ftp.ncbi.nih.gov/blast/db/

ascp（conda安装即可）高速下载链接：

1
ascp -v -k 1 -T -l 1000m -i /home/xxxxx/miniconda3/envs/bioinfTools/etc/asperaweb_id_dsa.openssh anonftp@ftp.ncbi.nih.gov:/blast/db/nt.00.tar.gz ./

需要将所有的文nt库文件全部下载！下载完成后解压即可。下载好的数据是已经完成建库的，可以直接运行blastn进行比对。

提取分类信息

处理blast结果

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
library(tidyverse)

vroom::vroom("./ncbi.nt.blast.result.ljy.txt", col_names = FALSE) -> df.ncbi.res

df.ncbi.res %>% 
  dplyr::filter(X3 == 100) -> df.ncbi.res.100 # 一共是258个

df.ncbi.res %>% 
  dplyr::filter(X3 != 100, X3 > 98)  %>% 
  dplyr::filter(!X1 %in% df.ncbi.res.100$X1) -> df.ncbi.res.over.98 # 一共是232个

df.ncbi.res.100 %>% 
  rbind(df.ncbi.res.over.98) -> df.ncbi.res.filtered

writexl::write_xlsx(df.ncbi.res.filtered, path = "./ncbi.nt.筛选后的结果.xlsx")

# 提取单独的NCBI accession
df.ncbi.res.filtered %>% 
  dplyr::select(2) %>% 
  dplyr::distinct_all()   %>% 
  magrittr::set_names(c("accession.version"))  %>% 
  data.table::as.data.table() -> id

id %>% 
  data.table::fwrite("./accession.id.txt")

准备accession2taxoid表格

从https://ftp.ncbi.nih.gov/pub/taxonomy/accession2taxid/下载核酸或者是蛋白的`accession2taxi`表格，如果是`gene bank的话需要下载nucl_gb.accession2taxid.gz`才行，整理成这个表格。

1
2
3
4
5
6
7
8
9
10
accession.version taxid
A00001.1 10641
A00002.1 9913
A00003.1 9913
A00004.1 32630
A00005.1 32630
A00006.1 32630
A00008.1 32630
A00009.1 32630
A00010.1 32630

将ncbi得到的结果与上述表格进行比较筛选

1
grep -f id.txt /home/xxx/database/nt/taxnonmy/nucl.gb.acc2taxid.txt > accession2taxoid.txt
1
head accession2taxoid.txt
1
2
3
4
5
6
7
8
9
10
AB001448.1 319
AB008001.1 303
AB013253.1 114707
AB015575.1 86473
AB019037.1 83814
AB022911.1 795097
AB023784.1 74315
AB024305.1 89568
AB038136.1 303
AB051692.1 143811

提取唯一的taxonomy编号

1
awk '{print $2}' accession2taxoid.txt | uniq -> taxo.id.txt

提取分类信息

1
taxonkit lineage taxo.id.txt |taxonkit reformat > taxo.id.2.taxonomy.txt

处理分类结果

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42

# 提取单独的NCBI accession
df.ncbi.res.filtered %>% 
  dplyr::select(2) %>% 
  dplyr::distinct_all()   %>% 
  magrittr::set_names(c("accession.version"))  %>% 
  data.table::as.data.table() -> id

id %>% 
  data.table::fwrite("./accession.id.txt")


# 读取taxonomy结果
vroom::vroom("./taxo.id.2.taxonomy.txt", col_names = FALSE)  %>% 
  dplyr::select(1,3)  %>% 
  dplyr::rename(taxoid = X1)  %>% 
  dplyr::mutate(界 = stringr::str_split(X3, ";") %>% sapply("[",1),
                门 = stringr::str_split(X3, ";") %>% sapply("[",2),
                纲 = stringr::str_split(X3, ";") %>% sapply("[",3),
                目 = stringr::str_split(X3, ";") %>% sapply("[",4),
                科 = stringr::str_split(X3, ";") %>% sapply("[",5),
                属 = stringr::str_split(X3, ";") %>% sapply("[",6),
                种 = stringr::str_split(X3, ";") %>% sapply("[",7))  %>% 
  dplyr::select(-2) -> df.taxo.res

# 读取accession id和taxo id表
vroom::vroom("./accession2taxoid.txt", col_names = FALSE)  %>% 
  magrittr::set_names(c("accession.id", "taxoid"))  %>% 
  dplyr::left_join(df.taxo.res, by = "taxoid")  %>% 
  dplyr::distinct_all() -> df.acc.taxo.res

# 合并blast结果和taxonomy结果
df.ncbi.res.filtered %>% 
  dplyr::select(1:3) %>% 
  magrittr::set_names(c("sample","accession.id","相似度")) %>% 
  dplyr::left_join(df.acc.taxo.res, by = "accession.id")  %>% 
  dplyr::mutate(temp = paste0(sample, 属))  %>% 
  dplyr::filter(!duplicated(temp))  %>% 
  dplyr::select(-temp) -> df.res.final

df.res.final %>% 
  writexl::write_xlsx("./最终结果.xlsx")
生物信息学
#生物信息学
NCBI-NT库本地化构建流程
https://lixiang117423.github.io/article/ncbint/
作者
小蓝哥
发布于
2023年2月22日
许可协议