はじめに

学生時代は主にPerl(Bioperl)を利用してバイオインフォマティクス分野の研究に携わっていましたが、最近はすっかりPythonが主流になっていますね。自分自身も業務でPythonを利用する機会が増えているので、Pythonの勉強がてらに学生時代を思い出しつつ、たまに「Python×バイオインフォマティクス」を学習してます。
今回はPython(Biopython)を利用した、Genbank形式ファイルをコーディング領域(CDS)のNucleotideまたはProteinのFASTAファイルへ変換する自作プログラムを紹介します。

プログラム入出力

入力Genbank形式ファイル
バイオインフォやってる人はお馴染みの形式。
テキストをそのまま貼り付けると長くなるので、内容についてはこちらから参照してください。
ファイル(sequence.gb)のダウンロードリンクはこちら (ファイルサイズは約10kb)
出力FASTA形式ファイル
上記のGenbank形式ファイルから下記のFASTA形式ファイルをそれぞれ出力可能としています。

Nucleotide配列FASTAファイル（sequence.fna）

>GENE00001 TCP1-beta
GATCCTCCATATACAACGGTATCTCCACCTCAGGTTTAGATCTCAACAACGGAACCATTGCCGACATGAGACAGTTAGGTATCGTCGAGAGTTACAAGCTAAAACGAGCAGTAGTCAGCTCTGCATCTGAAGCCGCTGAAGTTCTACTAAGGGTGGATAACATCATCCGTGCAAGACCAAGAACCGCCAATAGACAACATATGTAA
>GENE00002 Axl2p
ATGACACAGCTTCAGATTTCATTATTGCTGACAGCTACTATATCACTACTCCATCTAGTAGTGGCCACGCCCTATGAGGCATATCCTATCGGAAAACAATACCCCCCAGTGGCAAGAGTCAATGAATCGTTTACATTTCAAATTTCCAATGATACCTATAAATCGTCTGTAGACAAGACAGCTCAAATAACATACAATTGCTTCGACTTACCGAGCTGGCTTTCGTTTGACTCTAGTTCTAGAACGTTCTCAGGTGAACCTTCTTCTGACTTACTATCTGATGCGAACACCACGTTGTATTTCAATGTAATACTCGAGGGTACGGACTCTGCCGACAGCACGTCTTTGAACAATACATACCAATTTGTTGTTACAAACCGTCCATCCATCTCGCTATCGTCAGATTTCAATCTATTGGCGTTGTTAAAAAACTATGGTTATACTAACGGCAAAAACGCTCTGAAACTAGATCCTAATGAAGTCTTCAACGTGACTTTTGACCGTTCAATGTTCACTAACGAAGAATCCATTGTGTCGTATTACGGACGTTCTCAGTTGTATAATGCGCCGTTACCCAATTGGCTGTTCTTCGATTCTGGCGAGTTGAAGTTTACTGGGACGGCACCGGTGATAAACTCGGCGATTGCTCCAGAAACAAGCTACAGTTTTGTCATCATCGCTACAGACATTGAAGGATTTTCTGCCGTTGAGGTAGAATTCGAATTAGTCATCGGGGCTCACCAGTTAACTACCTCTATTCAAAATAGTTTGATAATCAACGTTACTGACACAGGTAACGTTTCATATGACTTACCTCTAAACTATGTTTATCTCGATGACGATCCTATTTCTTCTGATAAATTGGGTTCTATAAACTTATTGGATGCTCCAGACTGGGTGGCATTAGATAATGCTACCATTTCCGGGTCTGTCCCAGATGAATTACTCGGTAAGAACTCCAATCCTGCCAATTTTTCTGTGTCCATTTATGATACTTATGGTGATGTGATTTATTTCAACTTCGAAGTTGTCTCCACAACGGATTTGTTTGCCATTAGTTCTCTTCCCAATATTAACGCTACAAGGGGTGAATGGTTCTCCTACTATTTTTTGCCTTCTCAGTTTACAGACTACGTGAATACAAACGTTTCATTAGAGTTTACTAATTCAAGCCAAGACCATGACTGGGTGAAATTCCAATCATCTAATTTAACATTAGCTGGAGAAGTGCCCAAGAATTTCGACAAGCTTTCATTAGGTTTGAAAGCGAACCAAGGTTCACAATCTCAAGAGCTATATTTTAACATCATTGGCATGGATTCAAAGATAACTCACTCAAACCACAGTGCGAATGCAACGTCCACAAGAAGTTCTCACCACTCCACCTCAACAAGTTCTTACACATCTTCTACTTACACTGCAAAAATTTCTTCTACCTCCGCTGCTGCTACTTCTTCTGCTCCAGCAGCGCTGCCAGCAGCCAATAAAACTTCATCTCACAATAAAAAAGCAGTAGCAATTGCGTGCGGTGTTGCTATCCCATTAGGCGTTATCCTAGTAGCTCTCATTTGCTTCCTAATATTCTGGAGACGCAGAAGGGAAAATCCAGACGATGAAAACTTACCGCATGCTATTAGTGGACCTGATTTGAATAATCCTGCAAATAAACCAAATCAAGAAAACGCTACACCTTTGAACAACCCCTTTGATGATGATGCTTCCTCGTACGATGATACTTCAATAGCAAGAAGATTGGCTGCTTTGAACACTTTGAAATTGGATAACCACTCTGCCACTGAATCTGATATTTCCAGCGTGGATGAAAAGAGAGATTCTCTATCAGGTATGAATACATACAATGATCAGTTCCAATCCCAAAGTAAAGAAGAATTATTAGCAAAACCCCCAGTACAGCCTCCAGAGAGCCCGTTCTTTGACCCACAGAATAGGTCTTCTTCTGTGTATATGGATAGTGAACCAGCAGTAAATAAATCCTGGCGATATACTGGCAACCTGTCACCAGTCTCTGATATTGTCAGAGACAGTTACGGATCACAAAAAACTGTTGATACAGAAAAACTTTTCGATTTAGAAGCACCAGAGAAGGAAAAACGTACGTCAAGGGATGTCACTATGTCTTCACTGGACCCTTGGAACAGCAATATTAGCCCTTCTCCCGTAAGAAAATCAGTAACACCATCACCATATAACGTAACGAAGCATCGTAACCGCCACTTACAAAATATTCAAGACTCTCAAAGCGGTAAAAACGGAATCACTCCCACAACAATGTCAACTTCATCTTCTGACGATTTTGTTCCGGTTAAAGATGGTGAAAATTTTTGCTGGGTCCATAGCATGGAACCAGACAGAAGACCAAGTAAGAAAAGGTTAGTAGATTTTTCAAATAAGAGTAATGTCAATGTTGGTCAAGTTAAGGACATTCACGGACGCATCCCAGAAATGCTGTGA
>GENE00003 Rev7p
ATGAATAGATGGGTAGAGAAGTGGCTGAGGGTATACTTAAAATGCTACATTAATTTGATTTTATTTTATAGAAATGTATACCCACCTCAGTCATTCGACTACACTACTTACCAGTCATTCAACTTGCCGCAGTTCGTTCCCATTAATAGGCATCCTGCTTTAATTGACTATATAGAAGAACTTATACTGGATGTTCTTTCTAAATTAACGCACGTTTACAGATTTTCCATCTGCATTATTAATAAAAAGAACGATTTATGCATTGAAAAATACGTTTTAGATTTTAGTGAATTACAACATGTGGATAAAGACGATCAGATCATTACGGAAACTGAAGTGTTCGACGAATTCCGATCTTCCTTAAATAGTTTGATTATGCATTTGGAGAAATTACCTAAAGTCAACGATGACACAATAACATTTGAAGCAGTTATTAATGCGATCGAATTGGAACTAGGACATAAGTTGGACAGAAACAGGAGGGTCGATAGTTTGGAGGAAAAAGCAGAAATTGAAAGGGATTCAAACTGGGTTAAATGTCAAGAAGATGAAAATTTACCAGACAATAATGGTTTTCAACCTCCTAAAATAAAACTCACTTCTTTAGTCGGTTCTGACGTGGGGCCTTTGATTATTCATCAGTTTAGTGAAAAATTAATCAGCGGTGACGACAAAATTTTGAATGGAGTGTATTCTCAATATGAAGAGGGCGAGAGCATTTTTGGATCTTTGTTTTAA

Protein配列FASTAファイル（sequence.faa）

>GENE00001_AAA98665.1 TCP1-beta
SSIYNGISTSGLDLNNGTIADMRQLGIVESYKLKRAVVSSASEAAEVLLRVDNIIRARPRTANRQHM
>GENE00002_AAA98666.1 Axl2p
MTQLQISLLLTATISLLHLVVATPYEAYPIGKQYPPVARVNESFTFQISNDTYKSSVDKTAQITYNCFDLPSWLSFDSSSRTFSGEPSSDLLSDANTTLYFNVILEGTDSADSTSLNNTYQFVVTNRPSISLSSDFNLLALLKNYGYTNGKNALKLDPNEVFNVTFDRSMFTNEESIVSYYGRSQLYNAPLPNWLFFDSGELKFTGTAPVINSAIAPETSYSFVIIATDIEGFSAVEVEFELVIGAHQLTTSIQNSLIINVTDTGNVSYDLPLNYVYLDDDPISSDKLGSINLLDAPDWVALDNATISGSVPDELLGKNSNPANFSVSIYDTYGDVIYFNFEVVSTTDLFAISSLPNINATRGEWFSYYFLPSQFTDYVNTNVSLEFTNSSQDHDWVKFQSSNLTLAGEVPKNFDKLSLGLKANQGSQSQELYFNIIGMDSKITHSNHSANATSTRSSHHSTSTSSYTSSTYTAKISSTSAAATSSAPAALPAANKTSSHNKKAVAIACGVAIPLGVILVALICFLIFWRRRRENPDDENLPHAISGPDLNNPANKPNQENATPLNNPFDDDASSYDDTSIARRLAALNTLKLDNHSATESDISSVDEKRDSLSGMNTYNDQFQSQSKEELLAKPPVQPPESPFFDPQNRSSSVYMDSEPAVNKSWRYTGNLSPVSDIVRDSYGSQKTVDTEKLFDLEAPEKEKRTSRDVTMSSLDPWNSNISPSPVRKSVTPSPYNVTKHRNRHLQNIQDSQSGKNGITPTTMSTSSSDDFVPVKDGENFCWVHSMEPDRRPSKKRLVDFSNKSNVNVGQVKDIHGRIPEML
>GENE00003_AAA98667.1 Rev7p
MNRWVEKWLRVYLKCYINLILFYRNVYPPQSFDYTTYQSFNLPQFVPINRHPALIDYIEELILDVLSKLTHVYRFSICIINKKNDLCIEKYVLDFSELQHVDKDDQIITETEVFDEFRSSLNSLIMHLEKLPKVNDDTITFEAVINAIELELGHKLDRNRRVDSLEEKAEIERDSNWVKCQEDENLPDNNGFQPPKIKLTSLVGSDVGPLIIHQFSEKLISGDDKILNGVYSQYEEGESIFGSLF

実行コマンド・ソースコード

実行コマンド

実行ソースコードへ下記３つの引数を指定してプログラムを実行
① 入力Genbank形式ファイル、② 出力FASTA形式ファイル、③ 配列タイプ

出力コマンド例①：CDSのNucleotide配列FASTAファイル出力

$ python gbk2cds_fasta.py sequence.gb sequence.fna nucleotide

出力コマンド例②：CDSのProtein配列FASTAファイル出力

$ python gbk2cds_fasta.py sequence.gb sequence.faa protein

実行ソースコード

まずBiopythonのSeqIOでGenbank読み込み・構文解析を行い、CDSレコード情報のみを抽出しています。
その後、指定した配列タイプによりNucleotideまたはProtein配列を抜き出し、FASTA形式で書き込み実行しています。

gbk2cds_fasta.py

import sys
from typing import List

from Bio import SeqIO
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord


def main():
    # 引数設定
    args = sys.argv
    gbk_infile = args[1]
    cds_fasta_outfile = args[2]
    seqtype = args[3]

    # Genbank形式ファイル -> FASTA形式ファイル
    gbk2cds_fasta(gbk_infile, cds_fasta_outfile, seqtype)


def gbk2cds_fasta(
    gbk_infile: str, fasta_outfile: str, seqtype: str, id_prefix: str = "GENE"
) -> None:
    """GenBank形式ファイル内のCDS配列情報をFasta形式で出力

    Args:
        gbk_infile (str): 入力GenBank形式ファイル
        fasta_outfile (str): 出力Fasta形式ファイル
        seqtype (str): 出力配列タイプ("nucleotide" or "protein")
        id_prefix (str, optional): 出力配列連番IDのプレフィクス設定
    """
    gene_idx = 0
    cds_seq_record_list: List[SeqRecord] = []
    for record in SeqIO.parse(gbk_infile, "genbank"):
        # 翻訳配列レコードが存在するCDS配列を対象にFASTA出力
        cds_feature_list = [f for f in record.features if f.type == "CDS"]
        for cds_feature in cds_feature_list:
            qualifiers = cds_feature.qualifiers
            protein_id = qualifiers.get("protein_id", ["NA"])[0]
            product = qualifiers.get("product", ["NA"])[0]

            # 翻訳配列レコードが存在しない場合は出力スキップ
            if qualifiers.get("translation") is None:
                continue

            gene_idx += 1
            if seqtype == "nucleotide":
                seq_id = f"{id_prefix}{gene_idx:05d}"
                cds_seq = cds_feature.location.extract(record.seq)
            elif seqtype == "protein":
                seq_id = f"{id_prefix}{gene_idx:05d}_{protein_id}"
                cds_seq = Seq(qualifiers.get("translation")[0])
            else:
                raise KeyError(f"seqtype '{seqtype}' is not 'nucleotide|protein'")

            cds_seq_record = SeqRecord(seq=cds_seq, id=seq_id, description=product)
            cds_seq_record_list.append(cds_seq_record)

    # 配列行が1行のFASTA形式で出力
    SeqIO.write(cds_seq_record_list, fasta_outfile, "fasta-2line")


if __name__ == "__main__":
    main()

参考