Entrezと呼ばれるNCBIのデータベースからデータを取得するシステムをBiopythonを使って利用します。
やりたいこと
登録されている配列が解読された時に使用された、プライマーの配列情報をテキストファイルに書き出したいと思います。
簡単な流れ
- Biopythonを使って指定したIDからGenbank形式で情報を取得
- その情報から使用したプライマー情報を取得
- プライマー配列をGenbankIDのファイル名でテキストファイルに書き出す
注意点
情報を取得する際、サーバーに大きな負荷がかからないようにガイドラインに記載されているルールは守りましょう。
簡単にまとめると以下の通りです。2と3についてはBiopythonが勝手にやってくれるそうです
- 100以上のリクエストを送る場合は、週末かピークタイム以外の時間に行う
- 通常のアドレスではなく、http://eutils.ncbi.nlm.nih.gov/ を使用する
- リクエストの頻度を3リクエスト/1s 以下になるようにする
- e-mailとE-utilities APIの設定はしておく方がいい(Setting.txtなどにでも書き込んで実行する)
環境
- python3.8.5
- Ubuntu 20.04 LTS
- Biopythonはpip install biopythonなどを使ってダウンロードしておく
biopythonのモジュールなどをインポート
from Bio import Entrez, SeqIO
from Bio.SeqFeature import SeqFeature, FeatureLocation
from Bio.SeqRecord import SeqRecord
from pathlib import Path
import pathlib
NCBIに自己紹介をします。API keyはNCBIに登録して発行してもらっておきましょう。
# 自分のE-mailアドレスを記入
# Entrez.email = "A.N.Other@example.com"
# NCBIに登録したらAPI keyがもらえるのでそれを記入
# Entrez.api_key = "Your API key"
EFetch: Entrezからの完全な情報のダウンロード
NCBIのNucleotideデータベースのGenbankID: LC600705の情報を取得してみます。
NCBIで検索して情報を右上の’send to’からダウンロードするイメージで情報が取得できます。
# 取得したいデータのIDを指定
id = 'LC600705'
nucleotideのデータベースから、指定したIDの情報をGenbank形式で、テキストとして取得します。
handle = Entrez.efetch(db = 'nucleotide', id = id, rettype = 'gb', retmode = 'text')
print(handle.read())
LOCUS LC600705 808 bp DNA linear VRT 13-JAN-2021
DEFINITION Hemitrygon akajei CBM:ZF:20845 mitochondrial gene for 12S rRNA,
partial sequence.
ACCESSION LC600705
VERSION LC600705.1
KEYWORDS .
SOURCE mitochondrion Hemitrygon akajei (red stingray)
ORGANISM Hemitrygon akajei
Eukaryota; Metazoa; Chordata; Craniata; Vertebrata; Chondrichthyes;
Elasmobranchii; Batoidea; Myliobatiformes; Dasyatidae; Hemitrygon.
REFERENCE 1
AUTHORS Sado,T., Fukuchi,T. and Miya,M.
TITLE Reference data for MiFish metabarcoding analysis
JOURNAL Unpublished
REFERENCE 2 (bases 1 to 808)
AUTHORS Sado,T., Fukuchi,T. and Miya,M.
TITLE Direct Submission
JOURNAL Submitted (10-JAN-2021) Contact:Tetsuya Sado Natural History Museum
& Institute, Chiba; 955-2 Aoba-cho, Chuo-ku, Chiba, Chiba 260-8682,
Japan URL :http://www2.chiba-muse.or.jp/NATURAL/
FEATURES Location/Qualifiers
source 1..808
/organism="Hemitrygon akajei"
/organelle="mitochondrion"
/mol_type="genomic DNA"
/specimen_voucher="CBM:ZF:20845"
/db_xref="taxon:2704970"
/country="Japan:Chiba, Chiba, Chuo, Chiba Port Park, Beach
Plaza"
/PCR_primers="fwd_name: L-708-12S, fwd_seq:
ttayacatgcaagtatccgc, rev_name: H-1903-16S, rev_seq:
gtagctcgtytagtttcggg"
/note="synonym:Dasyatis akajei"
rRNA <1..>808
/product="12S ribosomal RNA"
ORIGIN
1 attccggtga gaacgcccta atcagcctac acatctaatt aggagctggt atcaggcaca
61 ctccaaagca gcccataaca cctcgctcgg ccacacccac aagggaattc agcagtgata
121 aacattgttc cataagcgca agcttgagtc aatcaaagtt ataaagagtt ggttaatctc
181 gtgccagcca ccgcggttat acgagtgacg caaactaata ttccacggcg ttaagggtga
241 ttagaaatat cttatccaaa ataaagttaa gaccccatta agctgtcata cgctctcatg
301 cttaaaaata tcactcacga aagtaacttt acataaaata gagtttttga cctcacgaca
361 gttaagaccc aaactaggat tagataccct actatgctta accgtaaaca tcgttataaa
421 taaatttacc ttaatacacc gcctgaacac tacaagcgct agcttaaaac ccaaaggact
481 tggcggtgct ccaaacccac ctagaggagc ctgttctata accgataatc cacgtttaac
541 ctcaccactt cttgcctttt tccgcctata taccgccgtc gtcagctcac cctgtgaaga
601 tacaacagta agcataatga ccttttaccc tcaacacgtc aggtcgaggt gtagcgaatg
661 aagtgggaag aaatgggcta cattctcttt ttagagtaca cgaacagaaa catgaaaaat
721 ttcttaaagg tggatttagc agtaagtaaa tttcagaata ttatactgaa accggctctg
781 gagcgcgcac acaccgcccg tcactctc
//
取得したGenbank形式のデータ情報をSeqIOを使って読み込みます。
次にSeqRecordを使ってパース(プログラムで扱えるようなデータ構造の集合体に変換)していきます。
例えば、要求したGenbankIDの情報の中から学名を表示するには以下のように記述します。
# Genbank形式で情報を取得
handle = Entrez.efetch(db = 'nucleotide', id = id, rettype = 'gb', retmode = 'text')
# 情報の読み込み
record = SeqIO.read(handle, "genbank")
# LC600705の学名を表示
organ = record.annotations['organism']
print(organ)
Hemitrygon akajei
Featuresの中の情報を取得するには以下のように記述する。
# Genbank形式で情報を取得
handle = Entrez.efetch(db = 'nucleotide', id = id, rettype = 'gb', retmode = 'text')
# 情報の読み込み
record = SeqIO.read(handle, "genbank")
# Feature情報を表示
feature = record.features
print(feature[0])
type: source
location: [0:808](+)
qualifiers:
Key: PCR_primers, Value: ['fwd_name: L-708-12S, fwd_seq: ttayacatgcaagtatccgc, rev_name: H-1903-16S, rev_seq: gtagctcgtytagtttcggg']
Key: country, Value: ['Japan:Chiba, Chiba, Chuo, Chiba Port Park, Beach Plaza']
Key: db_xref, Value: ['taxon:2704970']
Key: mol_type, Value: ['genomic DNA']
Key: note, Value: ['synonym:Dasyatis akajei']
Key: organelle, Value: ['mitochondrion']
Key: organism, Value: ['Hemitrygon akajei']
Key: specimen_voucher, Value: ['CBM:ZF:20845']
出力したい情報とそのキー(key)の対応が分かったので、欲しい情報を抜き出していきます。
# Genbank形式で情報を取得
handle = Entrez.efetch(db = 'nucleotide', id = id, rettype = 'gb', retmode = 'text')
# 情報の読み込み
record = SeqIO.read(handle, "genbank")
# typeがsorceの中でPCR_primerの情報を抜き出す
for feature in record.features:
if feature.type == 'source':
if 'PCR_primers' in feature.qualifiers:
primer_list = feature.qualifiers['PCR_primers']
# 文字列に変換
primer_str = str(primer_list[0])
# カンマで区切る
primer_split_list = primer_str.split(',')
# 各プライマー情報
FPrimer_name = primer_split_list[0].split(': ')[1]
FPrimer_seq = primer_split_list[1].split(': ')[1]
RPrimer_name = primer_split_list[2].split(': ')[1]
RPrimer_seq = primer_split_list[3].split(': ')[1]
# 出力
print(FPrimer_name)
print(FPrimer_seq)
print(RPrimer_name)
print(RPrimer_seq)
else:
primer_list = 'No information'
print(primer_list)
L-708-12S
ttayacatgcaagtatccgc
H-1903-16S
gtagctcgtytagtttcggg
取得したプライマーの情報をテキストファイルに書き出してみましょう。
書き出し先のフォルダを作成します。
# なければ作成 あれば作成しない
if not pathlib.Path('primer').exists():
pathlib.Path("primer").mkdir()
else:
print("primer folder exists")
先ほど作成したprimerフォルダの中にGenbankIDのファイル名でテキストファイルを作成してプライマー情報を書き出しします。
with open(f'./primer/{id}_primer.txt', mode = 'w') as primer:
primer.write(FPrimer_seq + '\n' + RPrimer_seq)
これでプライマーの情報が取得できるようになりました。
コードまとめ
少し改変したまとめコードになります。
作業ディレクトリにbiopython.pyというフォルダを作成し、下記コードをコピペして Ubuntu上で python3 biopython.py ________ (下線部にはGenbankID)と実行するれば、プライマーの情報がテキストファイルに出力されます。
プライマーの情報が無ければスクリプトが止まるようにしました。
emalとapi_keyは自身のものを記入してください。
# import module
from Bio import Entrez, SeqIO
from Bio.SeqFeature import SeqFeature, FeatureLocation
from Bio.SeqRecord import SeqRecord
from pathlib import Path
import pathlib
import sys
# Introduce myself
# Entrez.email = "A.N.Other@example.com"
# Entrez.api_key = "Your API key"
# Genbank ID
id = sys.argv[1]
# Get information by Genbank format.
handle = Entrez.efetch(db = 'nucleotide', id = id, rettype = 'gb', retmode = 'text')
record = SeqIO.read(handle, "genbank")
# Exract PCR_primer information from souce column
for feature in record.features:
if feature.type == 'source':
if 'PCR_primers' in feature.qualifiers:
primer_list = feature.qualifiers['PCR_primers']
# convert str
primer_str = str(primer_list[0])
# split comma
primer_split_list = primer_str.split(',')
# primer information
FPrimer_name = primer_split_list[0].split(': ')[1]
FPrimer_seq = primer_split_list[1].split(': ')[1]
RPrimer_name = primer_split_list[2].split(': ')[1]
RPrimer_seq = primer_split_list[3].split(': ')[1]
# print
print(FPrimer_name)
print(FPrimer_seq)
print(RPrimer_name)
print(RPrimer_seq)
else:
primer_list = 'Primer infromation does not exitsts.\nStop scripts.'
print(primer_list)
sys.exit(1)
# make folder
if pathlib.Path('primer').exists()==False:
pathlib.Path("primer").mkdir()
else:
print("primer folder exists")
# output primer folder
with open(f'./primer/{id}_primer.txt', mode = 'w') as primer:
primer.write(FPrimer_seq + '\n' + RPrimer_seq)
以上、ちょっとした趣味のお話でした。
参考
- Biopython Tutorial and Cookbook (http://biopython.org/DIST/docs/tutorial/Tutorial.html)
- Biopythonを利用したNCBIのEntrezデータベースへのアクセス(https://qiita.com/joemphilips/items/767c67524e4b7e328834)
- E-utilities API (https://ncbiinsights.ncbi.nlm.nih.gov/2017/11/02/new-api-keys-for-the-e-utilities/)