はじめに
学生時代は主に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()