|||
我们在做比较基因组学时会需要下载一些物种的同源基因序列,如果下载几十个菌株的序列,一个一个下载会很麻烦,所以这里我们进行批量下载,这里我们以解淀粉芽孢杆菌(bacillus amyloliquefaciens)为例。
首先在NCBI上搜索bacillus amyloliquefaciens:
在搜索结果里,我们可以选择查看该菌株的所有已测序菌株。
我们先看下该菌株的所有已测序的详细信息:在这个页面我们可以看到这些菌株的菌株名,样品名,接收号,以及很多基因组的统计信息。这里我们需要全基因组测序的菌株:
在点击搜索结果里面的 ftp下载链接 后,就进入到了如下图所示的ftp界面,我们需要从中找到下载链接的组成规律,
其实规律不难看出,就是由固定连接ftp://ftp.ncbi.nlm.nih.gov/genomes/genbank/bacteria/Bacillus_amyloliquefaciens/latest_assembly_versions/+菌株的Assembly号+菌株在NCBI上的编号+需要下载的文件名即可。
Assembly号我们可以从上面表格中找到,而这个编号在哪里呢,从下图我们可以知道,这个编号就在每个菌株的页面上,链接如下: https://www.ncbi.nlm.nih.gov/assembly/GCA_000196735.1
我们选择Compete后下载这个表格,然后使用命令行将strain列和Assembly列取出输出一个新的文件。
命令行如下:
cat genomes_proks848.txt |cut -f2,7 > bacillus_amyloiquefaciens_strain_list
## 显示刚刚生成的文件内容
less -S bacillus_amyloiquefaciens_strain_list
## 得到如下包含Strain和Assembly的列表
Strain Assembly
DSM 7 GCA_000196735.1
FS1092 GCA_004421045.1
Y2 GCA_000262385.1
KC41 GCA_008824205.1
SRCM101267 GCA_002173635.1
ZJU1 GCA_007362635.1
ALB69 GCA_003149695.1
SH-B74 GCA_003312895.1
ALB65 GCA_003149715.1
UMAF6639 GCA_001593765.1
ARP23 GCA_008245105.1
YP6 GCA_003868675.1
B15 GCA_001596755.1
UMAF6614 GCA_001593785.1
HK1 GCA_004006115.1
LL3 GCA_000204275.1
DH8030 GCA_008534375.1
LM2303 GCA_001889285.1
ALB79 GCA_003149795.1
MBE1283 GCA_001483885.1
Y14 GCA_001874385.1
H GCA_007827165.1
KHG19 GCA_000835145.1
X030 GCA_007923205.1
LFB112 GCA_000508265.1
XH7 GCA_000221645.1
TA208 GCA_000195515.1
ONU 553 GCA_008244925.1
WS-8 GCA_001922005.1
IT-45 GCA_000242855.2
S499 GCA_001586105.1
B-4 GCA_003390415.1
CC178 GCA_000494835.1
V417 GCA_008807995.1
V167 GCA_008808015.1
MT45 GCA_002209305.1
RD7-7 GCA_001705195.1
根据上面所说的,我们需要找到菌株在NCBI中的编号,这个编号是在网址: https://www.ncbi.nlm.nih.gov/assembly/GCA_000196735.1 中,该网址的源代码为(查看源代码方式时在网页上右击选择:查看网页源代码)
所以我们需要使用正则从网页中将编号提取出来,使用python代码实现如下:
import sys
import requests
import re
## 当输入的参数个数不等于2时,就打印本脚本的提示文本,以及参考输入文本,并且退出脚本
if len(sys.argv) != 2:
print('\nuseage:\npython3 get_assembly_name.py <input_txt> \nexample:\npython3 get_assembly_name.py bacillus_amyloiquefaciens_strain_list \n')
exit()
# 获取网页内容
def get_web_content(url):
header={'User-Agent':'Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/74.0.3729.169 Safari/537.36'}
response = requests.get(url,headers = header)
return response
# 打开bacillus_amyloiquefaciens_strain_list文件,取其中的Assembly和固定的网址组成组成新的网址,使用Requests包获取网页内容,然后使用正则获取我们想要的菌株编号
with open(sys.argv[1],'r') as read_file:
lines = read_file.readlines()
for line in lines[1:]:
col_list = line.strip().split('\t')
assembly_num = col_list[1].strip()
url = 'https://www.ncbi.nlm.nih.gov/assembly/' + assembly_num
web_content = get_web_content(url).text
assembly_name = re.search(r'<title>(.*)</title>',web_content).group(1)
assembly_name = assembly_name.split('-')[0].strip()
col_list.append(assembly_name)
print('\t'.join(col_list))
使用下面的命令得到包含菌株编号的新的文件
python get_assembly_name.py bacillus_amyloiquefaciens_strain_list > strain_list_add_assembly_name.txt
# 显示文件内容
less -S strain_list_add_assembly_name.txt
##
DSM_7 GCA_000196735.1 ASM19673v1
FS1092 GCA_004421045.1 ASM442104v1
Y2 GCA_000262385.1 ASM26238v1
KC41 GCA_008824205.1 ASM882420v1
SRCM101267 GCA_002173635.1 ASM217363v1
ZJU1 GCA_007362635.1 ASM736263v1
ALB69 GCA_003149695.1 ASM314969v1
SH-B74 GCA_003312895.1 ASM331289v1
ALB65 GCA_003149715.1 ASM314971v1
UMAF6639 GCA_001593765.1 ASM159376v1
ARP23 GCA_008245105.1 ASM824510v1
YP6 GCA_003868675.1 ASM386867v1
B15 GCA_001596755.1 ASM159675v1
UMAF6614 GCA_001593785.1 ASM159378v1
HK1 GCA_004006115.1 ASM400611v1
LL3 GCA_000204275.1 ASM20427v1
DH8030 GCA_008534375.1 ASM853437v1
LM2303 GCA_001889285.1 ASM188928v1
ALB79 GCA_003149795.1 ASM314979v1
MBE1283 GCA_001483885.1 ASM148388v1
Y14 GCA_001874385.1 ASM187438v1
H GCA_007827165.1 ASM782716v1
KHG19 GCA_000835145.1 ASM83514v1
X030 GCA_007923205.1 ASM792320v1
LFB112 GCA_000508265.1 ASM50826v1
XH7 GCA_000221645.1 ASM22164v1
TA208 GCA_000195515.1 ASM19551v1
ONU 553 GCA_008244925.1 ASM824492v1
WS-8 GCA_001922005.1 ASM192200v1
IT-45 GCA_000242855.2 ASM24285v2
S499 GCA_001586105.1 ASM158610v1
B-4 GCA_003390415.1 ASM339041v1
CC178 GCA_000494835.1 ASM49483v1
V417 GCA_008807995.1 ASM880799v1
V167 GCA_008808015.1 ASM880801v1
MT45 GCA_002209305.1 ASM220930v1
RD7-7 GCA_001705195.1 ASM170519v1
import sys
## 当输入的参数个数不等于3时,就打印本脚本的提示文本,以及参考输入文本,并且退出脚本
if len(sys.argv) != 3:
print('\nuseage:\npython3 download_ref_seq.py <strain_name> <input_txt> \nexample:\npython3 download_ref_seq.py Bacillus_amyloliquefaciens strain_list_add_assembly_name.txt \n')
exit()
## 打开strain_list_add_assembly_name.txt文件,按照下载的命令,取出文件的第二列Assembly号和第三列菌株编号进行组装,这里还增加了改名,在下载后就将文件名改为菌株名+原文件名的形式。
with open(sys.argv[2],'r') as read_file:
lines = read_file.readlines()
for line in lines:
line = line.strip().split('\t')
ID = line[1]+'_'+line[2]
file_name = ['_cds_from_genomic.fna.gz','_genomic.fna.gz','_genomic.gff.gz','_protein.faa.gz']
for i in file_name:
cmd_str = 'wget ftp://ftp.ncbi.nlm.nih.gov/genomes/genbank/bacteria/{0}/latest_assembly_versions/{1}/{2}'.format(sys.argv[1],ID,ID+i)
cmd_change_name = 'mv {0} {1}'.format(ID+i,line[0]+'_'+ID+i)
print(cmd_str)
print(cmd_change_name)
使用下面命令可以得到下载命令:
python download_ref_seq.py strain_list_add_assembly_name.txt > cmd_download_ref_seq.sh
得到的命令行文件内容如下:
接下来就是使用这个cmd_download_ref_seq.sh命令行文件进行下载序列文件了:
bash cmd_download_ref_seq.sh
最后得到的文件如下:
到这就完成了批量下载。
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-11-23 09:51
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社