Zenn
💻

Genbank形式のファイルからデータを取得して表形式にまとめる練習（Biopython）

に公開
Python
bioinformatics
biopython
tech

はじめに

Genbank形式のファイルからCDS feature内のqualifierの情報（のいくつか）を表形式のデータにまとめたい時があったのでBiopythonを使って処理するコードを検討してみます（主に原核生物が対象）。

動作確認 (Hardware)

-MacBook Pro
-チップ Apple M1
-macOS Ventura 13.6.6
-メモリ 16GB

動作確認 (Software)

-Python 3.13.3
-biopython 1.85
-pandas 2.2.3
-jupyter notebook

手順

1. コード例

各CDSの locus_tag、開始位置（start）、終了位置（end）、配列の向き（+/-）の情報（strand）、product、gene、pseudo の情報を取得するケースです。

'gene'のqualifierがない場合は'product'と同じ値を埋めるようにしています。

Python
#モジュールのインポート
from Bio import SeqIO
import pandas as pd

#GenBank形式のファイルを変数gb_fileに代入
gb_file = "/path_to/genbank_file"

#gb_fileをparseして、さらに辞書形式に変換してgb_recordに代入
gb_record = SeqIO.to_dict(SeqIO.parse(gb_file,"genbank"))

#空のデータフレームを用意
df_locus = pd.DataFrame([])
df_start = pd.DataFrame([])
df_end = pd.DataFrame([])
df_strand = pd.DataFrame([])
df_product = pd.DataFrame([])
df_gene = pd.DataFrame([])
df_pseudo = pd.DataFrame([])
df_CDS = pd.DataFrame([])

#データ取得
for i in gb_record:
    for f in gb_record[i].features:
        if f.type == "CDS":
            # qualifier keyに'/pseudogene'が含まれていればそのCDSの情報の取得をスキップ
            if 'pseudogene' in f.qualifiers:
                continue

            #locus_tagの名称を取得、qualifierがない場合は'--'で埋める
            locus = f.qualifiers.get('locus_tag', ['--'])[0]
            
            #CDSの開始位置を取得
            start = f.location.parts[0].start + 1
            
            #CDSの終了位置を取得
            end = f.location.parts[-1].end
            
            #+鎖か-鎖かの情報を取得(1 or -1)、1ならば +鎖、-1ならば -鎖
            strand = f.location.strand
            #1ならば "+"、-1ならば "-"に変換してstrand_strに代入
            if strand == 1:
                strand_str = "+"
            elif strand == -1:
                strand_str = "-"
            
            #productの情報を取得、qualifierがない場合は'--'で埋める
            product = f.qualifiers.get('product', ['--'])[0]
            
            #'/gene'の情報がなければproductと同じ値で埋める
            try:
                gene = f.qualifiers['gene']
            except KeyError:
                gene = product

            #'/pseudo'の情報の有無を検索
            if 'pseudo' in f.qualifiers:
                pseudo = str("pseudo")
            else:
                pseudo = str("--")

            #取得した値を各データフレームに追加            
            df_locus = pd.concat([df_locus, pd.DataFrame([locus])], ignore_index=True, axis=0)
            df_start = pd.concat([df_start, pd.DataFrame([start])], ignore_index=True, axis=0)
            df_end = pd.concat([df_end, pd.DataFrame([end])], ignore_index=True, axis=0)
            df_strand = pd.concat([df_strand, pd.DataFrame([strand_str])], ignore_index=True, axis=0)
            df_product = pd.concat([df_product, pd.DataFrame([product])], ignore_index=True, axis=0)
            df_gene = pd.concat([df_gene, pd.DataFrame([gene])], ignore_index=True, axis=0)
            df_pseudo = pd.concat([df_pseudo, pd.DataFrame([pseudo])], ignore_index=True, axis=0)

        #列方向に結合
        df_CDS =  pd.concat([df_locus, df_start, df_end, df_strand, df_product, df_gene, df_pseudo], axis=1)

#列名を設定
df_CDS.columns=["Locus", "Start", "End", "Strand", "Product", "Gene", "Pseudo"]

#全ての列の要素が一致している行があれば除く
df_CDS = df_CDS.drop_duplicates()

#インデックスを再指定
df_CDS.reset_index(inplace=True, drop=True)

2. 動作テスト

NCBIのデータを使っていくつかテストしてみます。

3. csvファイルに書き出す

csvファイルに保存する場合の例：

Python
df_CDS.to_csv('{}_annotation_table.csv'.format(gb_file), index=False)

カレントディレクトリに (genbank_file名)_annotation_table.csv ファイルができます。

おわりに

留意点：

  • 完全長になっていないドラフトゲノム（multi-contigの状態）にアノテーションをつけたファイルでは、コンティグが変わるごとにstartとendの値が1からにリセット
  • ひとつのCDSが複数のexonの結合のケースがある（真核生物での'join'表記など）、その場合startとendは不連続な領域の両端の値のみ取得

参考：

Discussion