从NCBI 上下载 gbff 文件并得到 CDS 信息

最近入门学习生信，面临的一个问题就是如何下载和预处理海量的基因组数据，尤其是 NCBI 上的 gbff 基因组注释文件，网络上翻了很久也找到对新手比较友好的代码，索性自己写了一段~

要处理的有两个问题，在此分别讨论：

一. 如何下载指定的文件（压缩包）且压缩包不损坏？

在此参考 NCBI 官方指南：Genomes Download FAQ (nih.gov)

个人之前是使用 wget 下载，常常出现压缩包损坏的情况，换了 rsync 没有再出现过该问题。

抓取文件的后缀（可以更换）：

- genomic.gbff.gz

- protein.faa.gz

网站地址（为啥我是用 https 而不是一步到位使用 rsync 或者 ftp？因为 https 我还可以直接进入该网页看看情况）：

https://ftp.ncbi.nlm.nih.gov/refseq/release/archaea/

https://ftp.ncbi.nlm.nih.gov/refseq/release/bacteria/

# 安装/导入 os、re、request 包
# 终端操作： pip install requests
import os
import re
import requests

def load_url_context(type ,url):
    # 从网页中提取指定文件目录
    request = requests.get(url)
    raw_list = re.compile(r'<a.*?>(.*?)</a>').finditer(request.text.strip())
    file_name = "_".join([type,"context.txt"])
    with open(file_name, "w") as f:
        for i in raw_list:
            x = i.group(1)
            if x.endswith("genomic.gbff.gz") or x.endswith("protein.faa.gz"):
                file_https = ''.join([url,x])
                # 使用 rsync 开头
                file_ftp = file_https.replace("https","rsync")
                #print(file_ftp)
                f.write(file_ftp)
                f.write('\n')
        f.close()    

if __name__ == "__main__":
    archaea_url = "https://ftp.ncbi.nlm.nih.gov/refseq/release/archaea/"
    bacteria_url = "https://ftp.ncbi.nlm.nih.gov/refseq/release/bacteria/"
    load_url_context("archaea", archaea_url)
    load_url_context("bacteria", bacteria_url)

上面的代码最后生成文件下载路径，再写一个小脚本开始下载数据 

#!/bin/bash

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

二. 是如何从注释文件中摘取自己想要的 CDS 注释信息

在上述网站上下载得到的数据是非常大的，尤其是细菌的基因组数据，但是有很多的注释信息我们用不到的，所以需要进一步操作精简一下内容。

在此我使用了 Biopython，用于 gbff 文件的导入，详情可以查看以下连接。

中文版：第5章 序列输入和输出 — Biopython-cn 0.1 文档

官网：Introduction to SeqIO · Biopython

本文我主要是想得到以下信息
 

Locus_tag	ORFStart..ORFStop	Strand	OrganismID	ContigID	Accession
KQY27_RS09175	167785:168877	-	Methanobrevibacter sp. TMH8 k121_54413-GCF_020148105.1	NZ_JAHLZE010000038.1	WP_224426279.1

import os
# 提前下载 Biopython 包
from Bio import SeqIO

# 搜索 gbff 文件，返回 gbff 文件路径
def GbffContext(path, types):
    file_path  = path + types
    files_path = os.path.abspath(file_path)
    all_files = os.listdir(files_path)
    context = []
    for file in all_files:
        if file.endswith("gbff"):
            gbff_file_path = file_path + '/' + file
            context.append(gbff_file_path)

    return context



#产生 cds.pty 文件
def CdsCoordinate(path, type):
    record_iterator = SeqIO.parse(path, 'genbank')
    cds_file = '_'.join([type,"cds.pty"])

    with open(cds_file,'a') as f:
        try:
            while(True):
                record = next(record_iterator)
                t = record.features
                for t_cds in t:
                    # cds 注释中必然含有 protein_id 的 key
                    if t_cds.qualifiers.__contains__('protein_id'):
                        f.write(str(t_cds.qualifiers["locus_tag"])[2:-2])# locus_tag, [2:-2]为了去掉括号和引号，下文同理
                        f.write('\t')

                        f.write(str(t_cds.location)[1:-4]) # ORFStart..ORFStop
                        f.write('\t')

                        f.write(str(t_cds.location)[-2])# Strand
                        f.write('\t')

                        f.write(record.description.split(',')[0])
                        f.write('-')
                        f.write(str(record.dbxrefs[-1])[9:]) # OrganismID
                        f.write('\t')

                        f.write(record.id)# ContigID
                        f.write('\t')
                        
                        f.write(str(t_cds.qualifiers["protein_id"])[2:-2])# Accession number
                        f.write("\n")
        except StopIteration as e:
            print("已经完成", path, "CDS内容的搜索")
    f.close()


# 输入参数
if __name__ == "__main__":
    path = 'your_path'
    type = 'archaea/ bacteria'

    gbff_context = GbffContext(path, type)
    for gbff in gbff_context:  
        print(gbff)
        CdsCoordinate(gbff, type)

总结

虽然代码效率很慢，but it works.继续苦逼的生信学习~