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などにでも書き込んで実行する)
参考
- 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/)
環境
- python3.8.5
- Ubuntu 20.04 LTS
- Biopythonはpip install biopythonなどを使ってダウンロードしておく
biopythonのモジュールなどをインポート
|
1 2 3 4 5 |
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に登録して発行してもらっておきましょう。
|
1 2 3 4 5 |
# 自分の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’からダウンロードするイメージで情報が取得できます。
|
1 2 |
# 取得したいデータのIDを指定 id = 'LC600705' |
nucleotideのデータベースから、指定したIDの情報をGenbank形式で、テキストとして取得します。
|
1 2 |
handle = Entrez.efetch(db = 'nucleotide', id = id, rettype = 'gb', retmode = 'text') print(handle.read()) |
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 |
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の情報の中から学名を表示するには以下のように記述します。
|
1 2 3 4 5 6 7 8 9 |
# Genbank形式で情報を取得 handle = Entrez.efetch(db = 'nucleotide', id = id, rettype = 'gb', retmode = 'text') # 情報の読み込み record = SeqIO.read(handle, "genbank") # LC600705の学名を表示 organ = record.annotations['organism'] print(organ) |
|
1 |
Hemitrygon akajei |
Featuresの中の情報を取得するには以下のように記述する。
|
1 2 3 4 5 6 7 8 9 |
# Genbank形式で情報を取得 handle = Entrez.efetch(db = 'nucleotide', id = id, rettype = 'gb', retmode = 'text') # 情報の読み込み record = SeqIO.read(handle, "genbank") # Feature情報を表示 feature = record.features print(feature[0]) |
|
1 2 3 4 5 6 7 8 9 10 11 |
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)の対応が分かったので、欲しい情報を抜き出していきます。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 |
# 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) |
|
1 2 3 4 |
L-708-12S ttayacatgcaagtatccgc H-1903-16S gtagctcgtytagtttcggg |
取得したプライマーの情報をテキストファイルに書き出してみましょう。
書き出し先のフォルダを作成します。
|
1 2 3 4 5 |
# なければ作成 あれば作成しない if not pathlib.Path('primer').exists(): pathlib.Path("primer").mkdir() else: print("primer folder exists") |
先ほど作成したprimerフォルダの中にGenbankIDのファイル名でテキストファイルを作成してプライマー情報を書き出しします。
|
1 2 |
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は自身のものを記入してください。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 |
# 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) |
以上、ちょっとした趣味のお話でした。@しばた
ABOUT ME
