🧬

biopythonのBio.SeqIOでgbffファイルを扱う方法

2022/03/16に公開約6,100字

はじめに

初投稿です。バイオインフォ系の記事はほとんど見かけないので場違い感がすごいですが、QiitaよりもUIがかっこいいのでZennで投稿させていただくことにしました。

この記事ではbiopythonのBio.SeqIOを使ってGenbankファイル(gbff)を操作する方法について、自分がよく使う部分だけを説明します。2022年に入ってから一か月おきくらいにこの要件が発生しては「どうやるんだっけ?」となってしまう自分のための解説記事です。

読込

基本的にgz圧縮されていると思うので、解凍込みの手順を示します。gzip -d等であらかじめ解凍しておいてももちろん構いません。

path = "./sample.gbff.gz"

import gzip
from Bio import SeqIO

with gzip.open(path, mode="rt", encoding="utf-8") as f:
    records = list(SeqIO.parse(f, "genbank"))
2022.03.17追記

イテレーターをジェネレーターと呼称して誤解を生みそうな文章だったので少し修正しました。
また、さらにSeqIO.parse()をいじっていて気づいたことですが、withでファイルがオープンされている間でなければnext()できないので、メモリ管理のために一つずつ読み込む場合には更に注意が必要そうです。

SeqRecordクラスオブジェクト

読み込んだrecordはSeqRecordクラスのインスタンスになります。

print(f"{records[0]=}")
print(f"{type(records[0])=}")

records[0]=SeqRecord(seq=Seq('GCGTTAAAACTTAAAGTGGAGGCGCGACCCGGAGTCGAACCGAGATCGACGGAT...ACA'), id='NZ_MCBT01000001.1', name='NZ_MCBT01000001', description='Shewanella colwelliana strain CSB03KR, whole genome shotgun sequence', dbxrefs=['BioProject:PRJNA224116', 'BioSample:SAMN05384437', 'Assembly:GCF_001735525.1'])
type(records[0])=<class 'Bio.SeqRecord.SeqRecord'>

読み込んだSeqRecordオブジェクトが個人的にとても取扱いにくい仕様に感じており、これ以降の処理が度々分からなくなります。

そのため、以降のセクションにて説明する内容がまとめの本質かもしれません。

インスタンス変数

以下のインスタンス変数はレコードの基本的な情報に参照する際に活用できます。強いて言えば最後のannotationsが重要で、それ以外はほとんど使用することはないかもしれません。

value 説明
dbxrefs アクセッション番号等 ['BioProject:PRJNA224116',(省略)]
id コンティグID 'NZ_MCBT01000001.1'
name シーケンスアクセッション 'NZ_MCBT01000001'
description コンティグの説明 'Shewanella colwelliana strain CSB03KR, whole genome shotgun sequence'
annotations 詳細なシーケンス情報 {'molecule_type': 'DNA', 'topology': 'linear',...(省略)

annotationsの下層

record.annotationにアクセスすると、非常に多くの情報がdict形式で得られます。

NCBIのエントリなので、しばしばメタ情報が他のサンプルと比較不能だったりするかもしれませんが、とりあえず今サンプルデータとして用いているものに含まれるメタ情報を羅列してみます。

value
molecule_type dna
topology linear
data_file_division CON
date 05-JUL-2021
accessions ['NZ_MCBT01000001', 'NZ_MCBT01000000']\
sequence_version 1
keywords ['WGS', 'RefSeq']
source Shewanella colwelliana
taxonomy ['Bacteria','Proteobacteria','Gammaproteobacteria','Alteromonadales','Shewanellaceae','Shewanella']
references 省略
comment 省略
contig 'join(MCBT01000001.1:1..325457)'
stductured_comment OrderedDictがネストしたオブジェクト

たとえば、taxonomyの情報が見たい時はdictのkeyで呼び出すようにやります。

records[0].annotations["taxonomy"]

SeqFeatureクラスオブジェクト

SeqFeatureクラスは遺伝子情報を記録したものになります。染色体1のstart, endの領域(strandは+/-)に、ほにゃららの機能を持つ遺伝子が存在している―――などといった情報です。

先ほどのSeqRecordクラスに対して、SeqRecord.featuresとやると、リスト形式でどさっとデータが出てきます。

たとえば、featureの中からrRNAだけを取り出したい場合にはこのようにします。

res = []
for record in records:
    minires = []
    for feature in record.features:
        if feature.type == "rRNA":
            minires.append(feature)
        else:
            pass
    if minires:
        res.append(minires)
    else:
        pass
res

[[SeqFeature(FeatureLocation(ExactPosition(149124), ExactPosition(150667), strand=1), type='rRNA')],
[SeqFeature(FeatureLocation(ExactPosition(241270), ExactPosition(241386), strand=-1), type='rRNA')],
[SeqFeature(FeatureLocation(ExactPosition(11), ExactPosition(1554), strand=1), type='rRNA'),
SeqFeature(FeatureLocation(ExactPosition(1909), ExactPosition(4803), strand=1), type='rRNA')],
[SeqFeature(FeatureLocation(ExactPosition(11), ExactPosition(1554), strand=1), type='rRNA'),
SeqFeature(FeatureLocation(ExactPosition(2403), ExactPosition(5297), strand=1), type='rRNA')],
[SeqFeature(FeatureLocation(BeforePosition(0), AfterPosition(516), strand=-1), type='rRNA')],
[SeqFeature(FeatureLocation(ExactPosition(319), AfterPosition(637), strand=-1), type='rRNA')],
[SeqFeature(FeatureLocation(ExactPosition(216), ExactPosition(332), strand=1), type='rRNA'),
SeqFeature(FeatureLocation(ExactPosition(497), ExactPosition(613), strand=1), type='rRNA')]]

SeqFeatureが保持する遺伝子の情報

ここで得られるSeqFeatureクラスのオブジェクトは一見単純な構造のデータ型に見えますが、print()してみると様々な情報を持っている事が分かります。

print(res[0][0])

type: rRNA
location: 149124:150667
qualifiers:
Key: db_xref, Value: ['RFAM:RF00177']
Key: inference, Value: ['COORDINATES: nucleotide motif:Rfam:12.0:RF00177', 'COORDINATES: profile:INFERNAL:1.1.1']
Key: locus_tag, Value: ['BEL05_RS03635']
Key: note, Value: ['Derived by automated computational analysis using gene prediction method: cmsearch.']
Key: old_locus_tag, Value: ['BEL05_09200']
Key: product, Value: ['16S ribosomal RNA']

sequenceの位置およびストランド(向き)は以下のようにアクセスできます。

print(res[0][0].location.start)
print(res[0][0].location.end)
print(res[0][0].location.strand)

149124
150667
1

qualifiersはdictとして情報が格納されています。
私は"product"にアクセスすることが多いです。その他motif情報などは活用できるか検討したことが無いので分かりません。

詳細は公式ドキュメントを見ていただいた方がよいかもしれません。

res[0][0].qualifiers['product']

['16S ribosomal RNA']

Seqクラスオブジェクト

seqオブジェクトは配列情報そのものですが、strクラスと同じような挙動になるように特殊メソッドを定義してくれてあります。

records[0].seq

Seq('GCGTTAAAACTTAAAGTGGAGGCGCGACCCGGAGTCGAACCGAGATCGACGGAT...ACA')

>>> len(records[0].seq)
325457
>>> str(records[0].seq)[:10]
'GCGTTAAAAC'
>>> records[0].seq[:10]
'GCGTTAAAAC'
>>> records[0].seq[:10] + "AAAAAAAAAAAAAAAA"
Seq('GCGTTAAAACAAAAAAAAAAAAAAAA')
>>> records[0].seq.count("AAG")
5692

挙動がstrクラスとほぼ一緒なのは良いのですが、この程度の操作だけならSeqクラスであるメリットが無く、私はすぐにstrクラスへ変換してしまいます。

多分アミノ酸翻訳などをやるとSeqクラスのありがたみがあるのだろうと思いますが、今のところそのような用事が無いのでテストしていません。

最後に

超シンプルな備忘録程度の情報しかありませんが、少なくとも未来の私は得すると思います。

ここまで読んでくれたあなたのためにもなれば嬉しく思います!

Discussion

ログインするとコメントできます