加载中...
根据assession number批量从NCB下载数据
发表于:2021-09-28 | 分类: 生物信息
字数统计: 1k | 阅读时长: 4分钟 | 阅读量:

有时候我们手里会得到一些 NCBI 的 assession number,且数量比较多,而我们真正需要的是序列,这时候手动挨个搜索和下载是不太现实的,除非是你闲得无事可做。其实有一个网页是可以批量下载序列的,即 https://www.ncbi.nlm.nih.gov/sites/batchentrez ,下面演示一下其用法。请就着文末视频食用。

  • 首先,准备一份列表文件,其中包含需要下载序列的 IDs,每行一个 ID。这里有一个从网上下载的 CaZY 数据库,本以为是序列文件,下载后才发现里面没有序列。这个文件包含三列,以制表符分隔各列,最后一列是 Assession number,因此前两列可以删掉。可以将文件内容复制到 Excel 中,删除前两列,将最后一列复制到一个新的文本文档中。也可以在支持正则表达式的文本编辑器中直接查找替换。刚刚的示例文件可以从这里下载。正则表达式查找的公式为 “.+\t (.+)”,其中 “.+” 代表的是任意多个字符,“\t” 匹配的是制表符,+ 是贪婪的,一直遇到最后一个 “\t” 才终止匹配。即 “.+\t” 匹配的是前两列以及第二列后面的制表符,最后的 “(.+)” 匹配的是第三列。小括号的作用是捕获匹配的内容。替换的公式为 “$1”,表示第一个小括号内的内容,即第三列。

  • 接下来将得到的列表文件提交至网站上以下载序列,需要选择对应的数据库,这里选择 protein,点击 “Retrieve” 开始下载。由于序列较多,因此反应比较慢,需要耐心等待。估计是崩了,再试一遍,文件包含 2650471 个 ID,估计服务器吃不消,实在不行就拆分成几个文件,分批次下载。我这里用的是 EmEditor 软件,按照 10000 行每个文件对整个文件进行了拆分,得到了 266 个文件,现在拿其中的一个做演示,看看服务器是否吃得消。看来一万个也太多,二十几个试了一下,莫得问题。方法就是酱紫的,至于可以一次性下载多少,各位自己去试吧。搞定!

  • 兄弟们不用试了,我已经试过了,一次只能搞几百个,对于几十万行的列表来说,手动逐个提交也是要命的,因此我写了一个 Perl 脚本 (download_NCBI.pl) 来完成该任务,不过只能在 Linux 下运行,代码如下:

#!/usr/bin/perl
use strict;
use warnings;
use LWP::Simple;
# Author: Liu hualin
# Date: Sep 28, 2021

# Usage: perl download_NCBI.pl 列表文件 序列类型(参照https://www.ncbi.nlm.nih.gov/sites/batchentrez数据库填写,常用的包括nucleotide, protein)

my @ids;
my $dbtype = $ARGV[1];# nucleotide, protein

system("split -l 300 $ARGV[0] splited_");

my @splited = glob("splited_*");
foreach  (@splited) {
	$_=~/splited_(.+)/;
	my $out = "seqs.$1.fasta";
	open IN, $_ || die;
	while (<IN>) {
		chomp;
		$_=~s/[\r\n]+//g;
		push @ids, $_;
	}
	getstore("http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=$dbtype&rettype=fasta&retmode=text&id=".join(",", @ids),"$out");
	@ids = ();
	close IN;
}

system("cat seqs.* > All.sequences.fas");
system("rm splited_* seqs.*");

运行方法也贼简单,将脚本和列表文件放在同一目录下,然后在 Linux 终端里输入如下命令:

perl download_NCBI.pl cazy_ids.txt protein

其中 cazy_ids.txt 为包含 assession number 的列表文件,protein 表示列表里的 ID 是蛋白。最后面的这个参数可以在 https://www.ncbi.nlm.nih.gov/sites/batchentrez 左上角的 Database 查询,但是要全部小写

运行一下,看看效果!

2000 years later...

2650471/300=8835 个文件,最终生成的序列文件名称为 “All.sequences.fas”,中间过程文件会被自动删除。千年以后拿结果,不管怎么说,总算解放了双手,让电脑慢慢去跑吧!

脚本获取

关注公众号 “生信之巅”,聊天窗口回复 “e7e9” 获取下载链接。

生信之巅微信公众号生信之巅小程序码

敬告:使用文中脚本请引用本文网址,请尊重本人的劳动成果,谢谢!Notice: When you use the scripts in this article, please cite the link of this webpage. Thank you!

上一篇:
Swissprot数据库的本地化与序列比对并与其他数据库快速mapping
下一篇:
Perl处理可恶的Windows换行符
本文目录
本文目录