🧬

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

2022/03/16に公開

Python

bio

tech

はじめに

初投稿です。バイオインフォ系の記事はほとんど見かけないので場違い感がすごいですが、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クラスのありがたみがあるのだろうと思いますが、今のところそのような用事が無いのでテストしていません。

最後に

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

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

はじめに

読込

SeqRecordクラスオブジェクト

インスタンス変数

annotationsの下層

SeqFeatureクラスオブジェクト

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

Seqクラスオブジェクト

最後に

Discussion

`SeqRecord`クラスオブジェクト

`annotations`の下層

`SeqFeature`クラスオブジェクト

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

`Seq`クラスオブジェクト