[Biopython] NCBI RefSeqの真核生物ゲノムから各遺伝子１本ずつタンパク質を取り出す

米国NCBI が公開するRefSeqは、DDBJ/ENA/GenBankといった公共配列データベースに登録された配列の中でも重要・高品質な参照配列 (Reference Sequences) を集めたもの。

ゲノム配列を収録したRefSeq genomesは、NCBIが独自にアノテーションを加え、遺伝子やその転写産物とタンパク質の情報を追加している。
NCBI RefSeqが収録する真核生物ゲノムは、NCBI Eukaryotic Genome Annotation Pipelineというパイプラインでアノテーションされている。

例：クルマエビ Penaeus japonicus
NCBI Penaeus japonicus Annotation Release 100

RefSeqのアノテーションに付属するタンパク質配列は、ゲノムアセンブリから欠落した配列がトランスクリプトーム等のデータで補われており、より信頼性が高い（と思われる）。

NCBI RefSeqの真核生物ゲノムアノテーション(GenBank形式ファイル*)から各遺伝子１本ずつ最長アイソフォームのタンパク質を抽出するスクリプトを書いた。

gb2nrfaa.py

Requirements

• Biopython

Usage

$ ./gb2nrfaa.py -h
usage: gb2nrfaa.py [-h] -i INPUT [-o OUTPUT] [-t TABLE]

Extract the longest isoforms of protein-coding genes from a NCBI RefSeq euaryotic genome assembly

optional arguments:
  -h, --help            show this help message and exit
  -i INPUT, --input INPUT
                        GenBank flat file format of the genomic sequence(s) (required)
  -o OUTPUT, --output OUTPUT
                        output fasta-formatted file (optional)
  -t TABLE, --table TABLE
                        output tab-separated table (optional)

Example:

$ wget https://ftp.ncbi.nlm.nih.gov/genomes/all/annotation_releases/27405/100/GCF_017312705.1_Mj_TUMSAT_v1.0/GCF_017312705.1_Mj_TUMSAT_v1.0_genomic.gbff.gz
$ gunzip GCF_017312705.1_Mj_TUMSAT_v1.0_genomic.gbff.gz
$ mv GCF_017312705.1_Mj_TUMSAT_v1.0_genomic.gbff Penaeus_japonicus.gb # optional
$ ./gb2nrfaa.py  -i Penaeus_japonicus.gb -o Penaeus_japonicus.faa # can take a few minutes
$ seqkit stats Penaeus_japonicus.faa
file                   format  type     num_seqs     sum_len  min_len  avg_len  max_len
Penaeus_japonicus.faa  FASTA   Protein    22,301  11,476,300       31    514.6   21,694

*FASTA+GFF形式ファイルだとpatchされた配列を反映できない。.faaファイルはLOCUS IDをヘッダに含んでおらず遺伝子情報が失われている。

2023-03-29追記: 各遺伝子とタンパク質配列の対応表をtab-seaprated fileとして出力できるようにした (default: out.tbl)。オプションは-tまたは--out-table。