搜索
查看: 2753|回复: 3

从NR库提取特定物种分类的序列

[复制链接]

29

主题

131

帖子

1200

积分

金牌会员

Rank: 6Rank: 6

积分
1200
发表于 2017-7-20 22:00:13 | 显示全部楼层 |阅读模式
本帖最后由 anlan 于 2017-7-20 22:11 编辑

1. 首先是软件的安装
下载NCBI的TaxonKit软件,http://bioinf.shenwei.me/taxonkit/download/,linux系统直接解压,接着
将taxonkit放到环境变量中
[AppleScript] 纯文本查看 复制代码
sudo cp taxonkit /usr/local/bin/

将两个文件(names.dmp和nodes.dmp)复制到用户目录下的隐藏文件夹.taxonkit中,
然后就能正常使用了
[AppleScript] 纯文本查看 复制代码
cp names.dmp ~/.taxonkit
cp nodes.dmp ~/.taxonkit


下载NCBI的csvtk软件,http://bioinf.shenwei.me/csvtk/download/,linux系统也是直接解压,即可使用

有一个数据文件需要下载,里面有NCBI的accession与taxid的对应关系,prot.accession2taxid.gz

2. 使用TaxonKit提取特定taxons下的所有taxid,以植物taxid 33090为例,--ids是指定taxid,--indent ""是将所列出的taxid左边的空格去除,以左对齐排列;然后也可以wc命令查看植物下有多少个taxid(有189153个taxids)

[AppleScript] 纯文本查看 复制代码
taxonkit list --ids 33090 --indent "" > plant.taxid.txt
wc -l plant.taxid.txt


3. 使用csvtk在prot.accession2taxid.gz文件中提取plant.taxid所有的accession,参数可以在http://bioinf.shenwei.me/csvtk/usage/查询到
[AppleScript] 纯文本查看 复制代码
zcat prot.accession2taxid.gz |csvtk -t grep -f taxid -P plant.taxid.txt |csvtk -t cut -f accession.version >plant.taxid.acc.txt


简单的说,感觉csvtk就是一个脚本集合工具,因此这步通过自己写的脚本也能实现的。大概1分钟后就能获得plant下所有的accession号,wc查下总共有8652634条

但是!!!如果只是想获得plant下的所有的accession号,其实上面3步都没必要做,因此NCBI网站上就可以直接下载,比如按照帖子http://www.biotrainee.com/thread-1818-1-1.html上的操作,在NCBI -> protein下检索txid33090[Organism] 即可显示出plant下所有的蛋白信息,然后单击Sent to -> File ->Accession List -> Creat File。所用时间可能还比上面的方法还要快。。。总共有8709559条,与上面的略有差异

但是你是其他大类物种,比如细菌等,并且还要提取好几个物种的accession号,那么上述的方法将会节约你不少时间

4. 如果是想建立NR子库,那么利用上面的得到的plant.taxid.acc.txt文件,使用blast+配套的程序即可,如:
[AppleScript] 纯文本查看 复制代码
blastdb_aliastool -gilist plant.taxid.acc.txt -db nr -out nr_plant -title nr_plant

5. 如果是想提取特定物种(比如植物)下的所有NR序列,那么可以按照http://bioinf.shenwei.me/taxonkit/tutorial/的方法,主要的也是利用blast+的blastdbcmd工具,如:
[AppleScript] 纯文本查看 复制代码
blastdbcmd -db nr -entry all -outfmt "%a\t%T" |csvtk -t grep -f 2 -P plant.taxid.acc.txt |csvtk -t cut -f 1 |blastdbcmd -db nr -entry_batch - -out nr.plant.fa

回复

使用道具 举报

1

主题

5

帖子

76

积分

注册会员

Rank: 2

积分
76
发表于 2017-7-27 05:52:41 | 显示全部楼层
本帖最后由 Sheng1982 于 2017-7-27 05:54 编辑

多谢楼主的好贴,对于我这种外行很有帮助。我想请教楼主,我下载了Nr到本地服务器,同时把植物的sequence.seq子库的accession列表下载了,用wc -l查看accession list有8711639条序列,但是用如下命令:../software/ncbi-blast-2.6.0+/bin/blastdb_aliastool -seqidlist sequence_plant.seq -db nr -out nr_plant -title nr_plant之后,只显示有5143766 sequences,这是正常情况吗,还是我这里出了问题?截图见附件。非常感谢。

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复 支持 反对

使用道具 举报

29

主题

131

帖子

1200

积分

金牌会员

Rank: 6Rank: 6

积分
1200
 楼主| 发表于 2017-7-28 19:51:00 | 显示全部楼层
Sheng1982 发表于 2017-7-27 05:52
多谢楼主的好贴,对于我这种外行很有帮助。我想请教楼主,我下载了Nr到本地服务器,同时把植物的sequence.s ...

我觉得的是正常的,因为从NCBI上下载的accession是NCBI上植物全部的accession,而NR提取出来的子库的是属于NR库的那部分植物的序列
回复 支持 反对

使用道具 举报

5

主题

26

帖子

286

积分

中级会员

Rank: 3Rank: 3

积分
286
发表于 2017-11-20 16:36:15 | 显示全部楼层
Sheng1982 发表于 2017-7-27 05:52
多谢楼主的好贴,对于我这种外行很有帮助。我想请教楼主,我下载了Nr到本地服务器,同时把植物的sequence.s ...

请问植物的sequence.seq子库的accession列表在哪里下载的?
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树    

GMT+8, 2019-3-24 03:28 , Processed in 0.039199 second(s), 24 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.