从NCBI refseq 中下载特定物种的蛋白质数据

前言

今天又和NCBI数据库干上了。由于NCBI奇奇怪怪的格式，导致我们下载特定物种（某个科、某个属）比较麻烦，手动一条条下载肯定是不现实的，而对于部分很少涉及干实验的生物研究人员来讲写代码也不容易。在此我分享一个思路，能够较为方便的下载数据：

工具

• Linux(wsl 也可以)
• python3
• rsync（NCBI推荐的下载方式）

代码与解析

我们使用bash来调用python脚本 RefseqDataDownload.py
RefseqDataDownload.py 需要两个物种分类学文件，数据来自于GTDB ,物种分类学文件可以自行准备或是使用附带文件A_Species_taxonomy.csv和B_Species_taxonomy.csv。请注意附带文件难以做到实时更新
RefseqDataDownload.py 需要三个参数：

1. 原核生物类型，输入 archaea/bacteria。
2. 下载数据的名称，例如属名Methanosarcina，注意如果是要下载某个种，请输入种属名(Methanosarcina_mazei，注意两者中间有下划线)。
3. 下载数据的分类学地位，其中 Domain=1, Phylum=2, Class=3, Order=4, Family=5, Genus=6, Species=7。

RefseqDataDownload.py 产生Contex.txt文件，包含对应物种的下载连接，进行接下来调用 rsync 下载

import os
import subprocess
import re
import argparse

class NameException(Exception): pass
subprocess.call("pip3 install pandas",shell=True)
subprocess.call("pip3 install numpy",shell=True)
subprocess.call("pip3 install requests",shell= True)

import pandas as pd
import numpy as np
import requests

def loadcsv(type):
    
    if type == 'archaea':
        csvname = 'A_Species_taxonomy.csv'
    elif type == 'bacteria':
        csvname = 'B_Species_taxonomy.csv'
    else:
        raise NameException("Error: please input archaea/bacteria!")
    
    SpeciesTaxonomy = pd.read_csv(csvname,names=['ID','Domain','Phylum','Class','Order','Family','Genus','Species'])
    return SpeciesTaxonomy

if __name__ == "__main__":
    ap = argparse.ArgumentParser(description = "按照分类地位下载NCBI Refseq的faa数据,数据来源例子:https://ftp.ncbi.nlm.nih.gov/genomes/refseq/archaea/Methanosarcina_mazei/representative/GCF_000007065.1_ASM706v1/GCF_000007065.1_ASM706v1_protein.faa.gz")
    ap.add_argument("-t", help = "原核生物的类型，输入 archaea/bacteria", required = True)
    ap.add_argument("-n", help = "下载数据的名称", required = True)
    ap.add_argument("-l", help = "分类等级, Domain=1, Phylum=2, Class=3, Order=4, Family=5, Genus=6, Species=7", required = True)

    opts = ap.parse_args()
    type = opts.t
    inputName = opts.n
    level = int(opts.l)

    SName = inputName.replace("_"," ")
    SpeciesTaxonomy = loadcsv(type)
    ColumnIndex = SpeciesTaxonomy.columns[level]
    IDinfo = (SpeciesTaxonomy.ID[SpeciesTaxonomy[ColumnIndex] == SName][0:][:3]== "GCF")

    OutputFrame = pd.DataFrame({'ID':SpeciesTaxonomy.ID[SpeciesTaxonomy[ColumnIndex] == SName],
    'Species':SpeciesTaxonomy.Species[SpeciesTaxonomy[ColumnIndex] == SName]})
    OutputFrame = OutputFrame[OutputFrame['ID'].apply(lambda x:x[0:3]) ==  'GCF']
    OutputFrame = OutputFrame.drop_duplicates(subset= ["Species"],keep='first')

    count = 0

    ContextFilename = inputName + 'Context.txt'
    with open(ContextFilename, "w") as f:
        for i in OutputFrame.Species:
            SpeciesName = i.replace(" ",'_')
            archaea_url = "https://ftp.ncbi.nlm.nih.gov/genomes/refseq/"+ type +"/"+SpeciesName+"/representative/"
            request = requests.get(archaea_url)
            raw_list = re.compile(r'<a.*?>(.*?)</a>').finditer(request.text.strip())

            for i in raw_list:
                x = i.group(1)
                if x[0:3] == 'GCF':
                    count = count + 1
                    faafile = x[:-1] + "_protein.faa.gz"
                    finalUrl_https = archaea_url + x + faafile
                    finalUrl_rsync = finalUrl_https.replace("https","rsync")
                    
                    f.write(finalUrl_rsync)
                    f.write('\n')
                    #print(finalUrl_rsync)

    print("完成搜索:",inputName,"\t序列总数:",count)
        

使用download.sh 完成下载：

#!/bin/bash

python3 RefseqDataDownload.py -t $1 -n $2 -l $3

#echo download data from NCBI

while read line
do 
    rsync --copy-links --recursive --times --verbose $line $2/
done < $2Context.txt

rsync下载速度一般，能翻墙的话速度会快些。

使用例子

# 在Linux系统上输入如下命令
bash download.sh archaea Methanosarcina 6
# 先产生一个MethanosarcinaContext.txt 文本文件
# 之后产生一个 Methanosarcina 文件夹。存放下载的 protein.faa.gz 数据

建议

如果想要下载其他的数据(.gbff，.fna，)可以修改RefseqDataDownload.py 代码，将faafile = x[:-1] + "_protein.faa.gz"改为对应的其他格式：fnafile = x[:-1] + "_genomic.fna.gz "
如果要下载原核生物的所有基因组数据，方法请参照 从NCBI 上下载 gbff 文件并得到 CDS 信息。

代码下载

https://download.youkuaiyun.com/download/LSD_1943/86711079?spm=1001.2014.3001.5501