Shiny从入门到入定——16-Escaping the graph

有时候我们手里会得到一些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!