Entrezによるデータベース検索システムを使って配列情報を取得します。 Entrezの検索システムは皆さんがNCBIで配列を検索したりするときに使うこのページそのものだと思っていただいたらいいと思います。
注意点
情報を取得する際、サーバーに大きな負荷がかからないようにガイドラインに記載されているルールは守る。
簡単にまとめると以下の通りです。2と3についてはBiopythonが勝手にやってくれるそうです。
- 100以上のリクエストを送る場合は、週末かピークタイム以外の時間に行う
- 通常のアドレスではなく、http://eutils.ncbi.nlm.nih.gov/ を使用する
- リクエストの頻度を3リクエスト/1s 以下になるようにする
- e-mailとE-utilities APIの設定はしておく方がいい (APIを登録しておくとリクエスト制限が緩和される) また、本記事はJupyter lab上でスクリプトを動かして、markdwon記法で記事を書いています。
Rstudioのような開発環境で、コードを実行すると結果がすぐわかるのはすごく使い勝手がよかったです。
動作環境
- Ubuntu 20.04 LTS
- python 3.8.5
- biopython 1.78
事前準備
Setting textファイルを作成
ディレクトリにSetting.txtを作成します。毎回パラメーターを記入するのは面倒なので、setting.txtに以下の情報をまとめて記入して、都度必要になったら呼び出すという形にします。
setting.txtの記載内容
emailアドレスとAPI keyは自身のものを記載してください。
####Caution####
All file names must not include space. Use "_" or "-" instead of space.
###############
# Setting values about all scripts.
Email = Your e-mail adress
Api_key = Your API key
# ==== Setting values about Download ID ===
# Variables for efetch
Database1 = nucleotide # Select the database to use.
Rettype1 = fasta # Select the format to download files.
Retmode1 = text # Select the file extension to download data.
Count_in_query1 = 20 # How many queries to request at a time. Setting range is 1 to 100000 at a time.
minday = 2000/01/01
Species = Plecoglossus altivelis
Genetic = mitochondrion
NCBIへの自己紹介
#emailアドレス
Email = "自分のメールアドレス"
#Apiキー
Api_key = "NCBIでアカウント登録した時に発行されたAPI番号"
emailやAPI keyは人によって異なる部分です。
パラメーター
添付画像の情報を指定していると思っていただけたらいいと思います。
#検索するデータベース
Database = "nucleotide" # ヌクレオチドのデータベースを使用
# データをダウンロードするときのフォーマット
Retmode = text # テキストで出力
RetmodeとRettypeの適切な組み合わせについてはこのページを参照ください。
# 1回の要求でダウンロードするデータ数
Count_in_query = 20
#いつ以降の配列を取得したいか
minday = 2000/01/01
# いつまでの配列を取得したいか
# 実行日を指定したい場合はスクリプト内に記載する為記載不要
maxday = 2021/01/27
# 対象種を指定
species = Plecoglossus altivelis
# 対象遺伝子を設定
genetic = mitochondrion
今回はアユ(Plecoglossus altivelis)のミトコンドリア遺伝子を指定します。
全体の設定
import module
import subprocess as sp
import pandas as pd
import pathlib, glob, re, pprint, os, datetime, shutil, time, sys
from Bio import Entrez, SeqIO
# from Bio.GenBank.Record import Record
from Bio.SeqFeature import SeqFeature, FeatureLocation
from Bio.SeqRecord import SeqRecord
from datetime import timedelta
from time import sleep
from dateutil.relativedelta import relativedelta
from subprocess import Popen, PIPE
from tqdm import tqdm
from termcolor import colored, cprint
結果をまとめるフォルダーのチェック
def check_results_folder():
# Reusltsフォルダの有無のチェック
if pathlib.Path("./Results").is_dir() == True:
shutil.rmtree("./Results") # フォルダを消して
pathlib.Path("./Results").mkdir(exist_ok=True) # 作成
else:
pathlib.Path("./Results").mkdir(exist_ok=True)
# IDをまとめるファイルのチェック
if pathlib.Path('./Results/ID.txt').is_file == True:
pathlib.Path("./Results/ID.txt").unlink()
else:
with open('./Results/ID.txt', mode='w') as i:
i.write("")
# Fastaファイルをまとめるファイルのチェック
if pathlib.Path("./Results/Fasta").is_dir() == False:
pathlib.Path("./Results/Fasta").mkdir(exist_ok=True)
else:
pass
# チェック
check_results_folder()
プログレスバーの色の設定
# Color setting
class Color:
BLACK = '\033[30m'
RED = '\033[31m'
GREEN = '\033[32m'
YELLOW = '\033[33m'
BLUE = '\033[34m'
PURPLE = '\033[35m'
CYAN = '\033[36m'
WHITE = '\033[37m'
END = '\033[0m'
BOLD = '\038[1m'
UNDERLINE = '\033[4m'
INVISIBLE = '\033[08m'
REVERCE = '\033[07m'
setting.txtの読み込み
with open('Setting.txt', encoding = "utf-8") as st:
txt = st.read()
NCBIへ自己紹介
reモジュールを使ってアドレスとAPI keyを取得します。
Entrez.email = re.search(r'Email\s*=\s*(\S+)', txt).group(1)
Entrez.api_key = re.search(r'Api_key\s*=\s*(\S+)', txt).group(1)
Genbank IDの取得
使用する変数の設定
# データベースを取得
db = re.search(r'Database1\s*=\s*(\S+)', txt).group(1)
# ファイルをダウンロードするときの形式
rettype = re.search(r'Rettype1\s*=\s*(\S+)', txt).group(1)
# データをダウンロードするときのフォーマット
retmode = re.search(r'Retmode1\s*=\s*(\S+)', txt).group(1)
# 1回の要求でダウンロードするデータ数
count_in_query = re.search(r'Count_in_query1\s*=\s*(\S+)', txt).group(1)
# いつからのデータを取得するか
minday = re.search(r'minday\s*=\s*(\S+)', txt).group(1)
# いつまでのデータを取得するか (常に今日になるようになっています)
maxday = datetime.date.today().strftime("%Y/%m/%d")
# 対象種を取得
species = re.search(r'Species\s*=\s*(.*)', txt).group(1)
# 対象遺伝子を設定
genetic= re.search(r'Genetic\s*=\s*(\S+)', txt).group(1)
esearchを動かす
esearchはEntrezの検索機能の部分で、efetchはダウンロード機能の部分です。IDを検索して取得したいのでesearchを使います。
# 配列情報の検索
handle_search = Entrez.esearch(db = db, retmax = count_in_query, term = f'{species}[Organism] AND {genetic}[filter]')
record_search = Entrez.read(handle_search)
# 取得したIDのList
SeqID = record_search["IdList"]
retmaxを20と指定したので、検索上位20本分の配列情報に付与されたGenbank IDを取得することができました。count_in_queryの数値は1~100000まで指定できるので、もっと大量の情報が欲しい時はこの数値を調整してみてください。
IDをテキストファイルに書き出す
# 操作しやすいようにID毎に改行する
str_ID = '\n'.join(SeqID)
# ResultsフォルダのID.txtに取得したIDを書き出す
with open('./Results/ID.txt', mode = 'a') as f:
f.write(str_ID + "\n")
これで、2000年1月1日以降に登録されたアユのミトコンドリアDNA配列に付与されているGenbankIDが取得できました。
次は取得したGenbankIDをもとにそれぞれの領域情報と、配列を取得していきます。
ResultsフォルダのID.txtには先ほど取得した20本分のGenbank IDが記載されています。それらを使ってそれぞれのIDの登録配列のDNA領域の情報を取得してみましょう。
# IDが記載されているID.txtの読み込み
ID_dataframe = pd.read_csv("./Results/ID.txt", encoding="UTF-8", header=None)
対象IDの配列を遺伝子領域ごとにFastaファイルに出力する
対象種はアユ、対象ゲノムはミトコンドリアDNAで検索した時にヒットした20個分のデータのGenbank IDがID.txtに格納されています。
ミトコンドリアDNAはリボソームRNAやtransfer RNAなど数十種類の遺伝子を有しています。MiFishプライマーなどは12s rRNA遺伝子に設計されていたり、系統解析に使われるようなシトクロムb遺伝子はミトコンドリアDNAにコードされています。
# New function define
def makeSeqRecord(seq_obj, id_for_fasta):
record = SeqRecord(seq_obj, id=id_for_fasta, description="")
return record
# rRNA or CDS or tRNA misc_feature
feature_str = 'rRNA'
# Genbank ID数を取得
id_num = len(ID_dataframe)
# Genbank ID の数だけForで情報を取得する
for i in tqdm(range(0, id_num), desc = Color.CYAN + 'Progress rate' + Color.END):
# IDの配列情報をGenbank形式で取得
handle_info = Entrez.efetch(db = db, id = ID_dataframe.iloc[i,0], rettype = 'gb', retmode = 'text')
for record_info in SeqIO.parse(handle_info, "genbank"): # gbを指定するとGenbank形式の情報を取得する
# リスト格納
record_list = []
# 配列全体を格納
parent_sequence = record_info.seq
# 学名を格納
organism = record_info.annotations['organism']
# Accession numberを格納
acc = record_info.id
# 各領域ごとに配列を切り分けて配分
for feature in record_info.features:
# 不明な遺伝子領域 misc_feature 以外の場合
if feature.type != 'misc_feature' and feature_str == feature.type:
sequence_slice = feature.location.extract(parent_sequence)
# 配列の末端取得
strand_int = feature.location.strand
if strand_int == 1:
strand = 'top' # 5'末端
elif strand_int == -1:
strand = 'btm' # 3'末端
else:
strand = 'cannotbe'
# 配列の開始位置指定
start = feature.location.parts[0].start.position + 1
# 配列の終了位置指定
end = feature.location.parts[-1].end.position
# 遺伝子領域の情報取得
if 'gene' in feature.qualifiers:
gene = feature.qualifiers['gene'][0]
else:
gene = 'NA'
# 転写産物？の情報取得
if 'product' in feature.qualifiers:
product = feature.qualifiers['product'][0]
else:
product = 'NA'
# Fasta形式にした時の情報部分(> ~~~ の部分)の作成
id_line = acc + '_' + organism + '_' + gene + '_' + product + '_'
id_line += ' (' + str(start) + '..' + str(end) + ', strand=' + strand + ')'
# レコードから取得した情報をFasta形式にする
record_list = makeSeqRecord(sequence_slice, id_line)
# Fastaファイルの出力
if feature_str != 'D-loop':
fasta_file_name = './Results/Fasta/' + acc + f'_{product}_' + feature.type + '.fasta'
else:
fasta_file_name = './Results/Fasta/' + acc + f'_{feature_str}_test.fasta'
# 一応0以上かチェック
type_count = len(record_list)
with open(fasta_file_name, mode = 'w') as ffn:
if type_count > 0:
SeqIO.write(record_list, ffn, 'fasta')
else :
ffn.write('No Annotation Found')
# merge fasta file name
mfa = f'./Results/Fasta/mergefa.fas'
with open(mfa, mode = 'a') as mfa:
if type_count > 0:
SeqIO.write(record_list, mfa, 'fasta')
else:
pass
with open("./Results/list.txt", mode = 'w') as spl:
spl.write(organism + '\n')
# 不明な遺伝子領域 misc_feature の場合(D-loopがここに含まれる場合があるため作成)
elif feature_str == 'misc_feature' and feature_str == feature.type or feature_str == 'misic_feature' and re.search('(d|D)-(l|L)oop', feature.type) == True:
# feature.qualifiers フィールドにD-loopの記載がある場合があるため
if 'note' in feature.qualifiers:
note = feature.qualifiers['note'][0]
else:
note = 'No Annotation'
sequence_slice = feature.location.extract(parent_sequence)
# 配列の末端
strand_int = feature.location.strand
if strand_int == 1:
strand = 'top' # 5'末端
elif strand_int == -1:
strand = 'btm' # 3'末端
else:
strand = 'cannotbe'
# 配列の開始位置
start = feature.location.parts[0].start.position + 1
# 配列の終了位置
end = feature.location.parts[-1].end.position
# 転写産物？の情報
product = 'NA'
# Fasta形式にした時の情報部分の作成
id_line = acc + '_' + organism + '_' + note + '_'
id_line += ' (' + str(start) + '..' + str(end) + ', strand=' + strand + ')'
# レコードから取得した情報をFasta形式にする
record_list = makeSeqRecord(sequence_slice, id_line)
# Fastaファイルの出力
fasta_file_name = './Results/Fasta/' + acc +'_' + feature.type + '.fasta'
type_count = len(record_list)
with open(fasta_file_name, mode = 'w') as ffn:
if type_count > 0:
SeqIO.write(record_list, ffn, 'fasta')
else :
ffn.write('No Annotation Found')
# merge fasta file name
mfa = './Results/Fasta/mergefa.fas'
with open(mfa, mode = 'a') as mfa:
if type_count > 0:
SeqIO.write(record_list, mfa, 'fasta')
else:
pass
# リストに追加された種名を記入
add_organism_list = []
add_organism_list = add_organism_list.append(organism)
print(add_organism_list)
with open("./Results/list.txt", mode = 'w') as spl:
spl.write(organism + '\n')
else:
pass
以下のフォルダ構成になります。1本しかhitしてないのでmergefa.fasには1本しかないですが、対象種や対象feature、retmaxを変えたりしたら複数配列が取得できます。
Fasta
- LC600708.1_12S ribosomal RNA_rRNA.fasta
- mergefa.fas
ID.txt
list.txt
DNA領域に構わず、対象IDの配列を一括取得する
領域ごとに切り出すとかファイルを分けるとかせずに、IDに対応する配列を全部取得するんだ。という方はこちらを使用してみてください。2-1の部分をごっそり以下のコードに変更します。
# New function define
def makeSeqRecord(seq_obj, id_for_fasta):
record = SeqRecord(seq_obj, id=id_for_fasta, description="")
return record
# IDが記載されているID.txtの読み込み
ID_dataframe = pd.read_csv("./Results/ID.txt", encoding="UTF-8", header=None)
# Genbank ID数を取得
id_num = len(ID_dataframe)
# GWnbank ID の数だけForで情報を取得する
for i in tqdm(range(0, id_num), desc = Color.CYAN + 'Progress rate' + Color.END):
handle = Entrez.efetch(db = db, id = ID_dataframe.iloc[i,0], rettype = 'gb', retmode = 'text')
for record_info in SeqIO.parse(handle, "genbank"): # gbを指定するとGenbank形式の情報を取得する
# リスト格納
record_list = []
# 配列全体を格納
parent_sequence = record_info.seq
# 学名を格納
organism = record_info.annotations['organism']
# Accession numberを格納
acc = record_info.id
# 各領域ごとに配列を切り分けて配分
for feature in record_info.features:
# 不明な遺伝子領域 misc_feature 以外の場合
if feature.type == 'source':
sequence_slice = feature.location.extract(parent_sequence)
# 配列の末端取得
strand_int = feature.location.strand
if strand_int == 1:
strand = 'top' # 5'末端
elif strand_int == -1:
strand = 'btm' # 3'末端
else:
strand = 'cannotbe'
# 配列の開始位置指定
start = feature.location.parts[0].start.position + 1
# 配列の終了位置指定
end = feature.location.parts[-1].end.position
# Fasta形式にした時の情報部分(> ~~~ の部分)の作成
id_line = acc + '_' + organism + '(' + str(start) + '..' + str(end) + ', strand=' + strand + ')'
# レコードから取得した情報をFasta形式にする
record_list = makeSeqRecord(parent_sequence, id_line)
# アクセッションナンバー.fa のファスタファイルでファイル名を指定
fasta_file_name = './Results/Fasta/' + record_info.id + '.fasta'
# 一応0以上かチェック
type_count = len(record_list)
# 書き出し
with open(fasta_file_name, mode = 'w') as ffn:
if type_count > 0:
SeqIO.write(record_list, ffn, 'fasta')
else :
ffn.write('No Annotation Found')
# merge fasta file name
mfa = './Results/Fasta/mergefa.fas'
with open(mfa, mode = 'a') as mfa:
if type_count > 0:
SeqIO.write(record_list, mfa, 'fasta')
else:
pass
以下がResultsフォルダ内の構成になります。mergefa.fasには20本分のfastaファイルが一つにまとめられています。
Fasta
- Fastaファイル20本
- mergefa.fas
ID.txt 場合に応じて使い分けてみてください。
番外編
2-2で取得した配列の特徴をmatplotlibで視覚化してみる。
# 成型
with open('./Results/Fasta/mergefa.fas', mode = 'r') as fa:
if i == 0:
fasta = fa.read()
fasta_rp = re.sub(r'>', r'\n>', fasta)
fasta_rp = re.sub(r"top\)\n", r"top)\t", fasta_rp)
fasta_rp = re.sub(r"btm\)\n", r"btm)\t", fasta_rp)
fasta_rp = re.sub(r'([A-Z]|-$)\n', r'\1', fasta_rp)
with open('./Results/Fasta/mergefa_edit.fas', mode = "w") as wr:
wr.write(fasta_rp)
else:
fasta = fa.read()
fasta_rp = re.sub(r'>', r'\n>', fasta)
fasta_rp = re.sub(r'^\n', '', fasta_rp)
fasta_rp = re.sub(r"top\)\n", r"top)\t", fasta_rp)
fasta_rp = re.sub(r"btm\)\n", r"btm)\t", fasta_rp)
fasta_rp = re.sub(r'([A-Z]|-$)\n', r'\1', fasta_rp)
# 書き出し
with open('./Results/Fasta/merge_edit.fas', mode='w') as wr:
wr.write(fasta_rp)
# テーブルデータとして読み込み
mfa = pd.read_csv('./Results/Fasta/merge_edit.fas', encoding="UTF-8", header = None, index_col = False, sep = "\t")
# 列名を変更
mfa = mfa.rename(columns={0:"ID",1:"seq"})
# 配列の長さを算出
mfa['seq_len'] = mfa['seq'].str.len()
# 配列長でhistグラムを描画してみる
import matplotlib.pyplot as plt
from pylab import rcParams
rcParams['figure.figsize'] = 20, 5
rcParams['font.size'] = 25
mfa.hist(bins=50)
plt.tight_layout()
plt.show()
まとめのコードなど
いくつかのパターンを同時並行で解説していったので、流れが分かりにくくなっていたと思います。
パターンに分けてコードをまとめて記載するので、その通りやれば類似した結果が得られると思います。
Setting.txtの内容
####Caution####
All file names must not include space. Use "_" or "-" instead of space.
###############
# Setting values about all scripts.
Email = 'Your e-mail adress'
Api_key = 'Your API key'
# ==== Setting values about Download ID ===
# Variables for efetch
Database1 = nucleotide # Select the database to use.
Rettype1 = fasta # Select the format to download files.
Retmode1 = text # Select the file extension to download data.
Count_in_query1 = 20 # How many queries to request at a time. Setting range is 1 to 100000 at a time.
minday = 2000/01/01
Species = Plecoglossus altivelis
Genetic = mitochondrion
Feature = rRNA
対象IDの配列を遺伝子領域ごとにFastaファイルに出力するスクリプト
# --- import module --- #
import subprocess as sp
import pandas as pd
import pathlib, glob, re, pprint, os, datetime, shutil, time, sys
from Bio import Entrez, SeqIO
# from Bio.GenBank.Record import Record
from Bio.SeqFeature import SeqFeature, FeatureLocation
from Bio.SeqRecord import SeqRecord
from datetime import timedelta
from time import sleep
from dateutil.relativedelta import relativedelta
from subprocess import Popen, PIPE
from tqdm import tqdm
from termcolor import colored, cprint
# --- フォルダーのチェック --- #
def check_results_folder():
# Reusltsフォルダの有無のチェック
if pathlib.Path("./Results").is_dir() == True:
shutil.rmtree("./Results") # フォルダを消して
pathlib.Path("./Results").mkdir(exist_ok=True) # 作成
else:
pathlib.Path("./Results").mkdir(exist_ok=True)
# IDをまとめるファイルのチェック
if pathlib.Path('./Results/ID.txt').is_file == True:
pathlib.Path("./Results/ID.txt").unlink()
else:
with open('./Results/ID.txt', mode='w') as i:
i.write("")
# Fastaファイルをまとめるファイルのチェック
if pathlib.Path("./Results/Fasta").is_dir() == False:
pathlib.Path("./Results/Fasta").mkdir(exist_ok=True)
else:
pass
# --- チェック --- #
check_results_folder()
# --- Color setting --- #
class Color:
BLACK = '\033[30m'
RED = '\033[31m'
GREEN = '\033[32m'
YELLOW = '\033[33m'
BLUE = '\033[34m'
PURPLE = '\033[35m'
CYAN = '\033[36m'
WHITE = '\033[37m'
END = '\033[0m'
BOLD = '\038[1m'
UNDERLINE = '\033[4m'
INVISIBLE = '\033[08m'
REVERCE = '\033[07m'
# --- setting.txtの読み込み --- #
with open('Setting.txt', encoding = "utf-8") as st:
txt = st.read()
# --- NCBIへの自己紹介 --- #
Entrez.email = re.search(r'Email\s*=\s*(\S+)', txt).group(1)
Entrez.api_key = re.search(r'Api_key\s*=\s*(\S+)', txt).group(1)
# --- setting.txtから変数の情報を取得する --- #
# データベースを取得
db = re.search(r'Database1\s*=\s*(\S+)', txt).group(1)
# ファイルをダウンロードするときの形式
rettype = re.search(r'Rettype1\s*=\s*(\S+)', txt).group(1)
# データをダウンロードするときのフォーマット
retmode = re.search(r'Retmode1\s*=\s*(\S+)', txt).group(1)
# 1回の要求でダウンロードするデータ数
count_in_query = re.search(r'Count_in_query1\s*=\s*(\S+)', txt).group(1)
# いつからのデータを取得するか
minday = re.search(r'minday\s*=\s*(\S+)', txt).group(1)
# いつまでのデータを取得するか (常に今日になるようになっています)
maxday = datetime.date.today().strftime("%Y/%m/%d")
# 対象種を取得
species = re.search(r'Species\s*=\s*(.*)', txt).group(1)
# 対象遺伝子を設定
genetic= re.search(r'Genetic\s*=\s*(\S+)', txt).group(1)
# --- IDを検索する --- #
handle_search = Entrez.esearch(db = db, retmax = count_in_query, term = f'{species}[Organism] AND {genetic}[filter]')
record_search = Entrez.read(handle_search)
# 取得したIDのList
SeqID = record_search["IdList"]
# 操作しやすいようにID毎に改行する
str_ID = '\n'.join(SeqID)
# ResultsフォルダのID.txtに取得したIDを書き出す
with open('./Results/ID.txt', mode = 'a') as f:
f.write(str_ID + "\n")
# Dataframe の形でIDの読み込み
ID_dataframe = pd.read_csv("./Results/ID.txt", encoding="UTF-8", header=None)
# New function define
def makeSeqRecord(seq_obj, id_for_fasta):
record = SeqRecord(seq_obj, id=id_for_fasta, description="")
return record
# rRNA or CDS or tRNA misc_feature
feature_str= re.search(r'Feature\s*=\s*(\S+)', txt).group(1)
# Genbank ID数を取得
id_num = len(ID_dataframe)
# Genbank ID の数だけForで情報を取得する
for i in tqdm(range(0, id_num), desc = Color.CYAN + 'Progress rate' + Color.END):
# IDの配列情報をGenbank形式で取得
handle_info = Entrez.efetch(db = db, id = ID_dataframe.iloc[i,0], rettype = 'gb', retmode = 'text')
for record_info in SeqIO.parse(handle_info, "genbank"): # gbを指定するとGenbank形式の情報を取得する
# リスト格納
record_list = []
# 配列全体を格納
parent_sequence = record_info.seq
# 学名を格納
organism = record_info.annotations['organism']
# Accession numberを格納
acc = record_info.id
# 各領域ごとに配列を切り分けて配分
for feature in record_info.features:
# 不明な遺伝子領域 misc_feature 以外の場合
if feature.type != 'misc_feature' and feature_str == feature.type:
sequence_slice = feature.location.extract(parent_sequence)
# 配列の末端取得
strand_int = feature.location.strand
if strand_int == 1:
strand = 'top' # 5'末端
elif strand_int == -1:
strand = 'btm' # 3'末端
else:
strand = 'cannotbe'
# 配列の開始位置指定
start = feature.location.parts[0].start.position + 1
# 配列の終了位置指定
end = feature.location.parts[-1].end.position
# 遺伝子領域の情報取得
if 'gene' in feature.qualifiers:
gene = feature.qualifiers['gene'][0]
else:
gene = 'NA'
# 転写産物？の情報取得
if 'product' in feature.qualifiers:
product = feature.qualifiers['product'][0]
else:
product = 'NA'
# Fasta形式にした時の情報部分(> ~~~ の部分)の作成
id_line = acc + '_' + organism + '_' + gene + '_' + product + '_'
id_line += ' (' + str(start) + '..' + str(end) + ', strand=' + strand + ')'
# レコードから取得した情報をFasta形式にする
record_list = makeSeqRecord(sequence_slice, id_line)
# Fastaファイルの出力
if feature_str != 'D-loop':
fasta_file_name = './Results/Fasta/' + acc + f'_{product}_' + feature.type + '.fasta'
else:
fasta_file_name = './Results/Fasta/' + acc + f'_{feature_str}_test.fasta'
# 一応0以上かチェック
type_count = len(record_list)
with open(fasta_file_name, mode = 'w') as ffn:
if type_count > 0:
SeqIO.write(record_list, ffn, 'fasta')
else :
ffn.write('No Annotation Found')
# merge fasta file name
mfa = f'./Results/Fasta/mergefa.fas'
with open(mfa, mode = 'a') as mfa:
if type_count > 0:
SeqIO.write(record_list, mfa, 'fasta')
else:
pass
with open("./Results/list.txt", mode = 'w') as spl:
spl.write(organism + '\n')
# 不明な遺伝子領域 misc_feature の場合(D-loopがここに含まれる場合があるため作成)
elif feature_str == 'misc_feature' and feature_str == feature.type or feature_str == 'misic_feature' and re.search('(d|D)-(l|L)oop', feature.type) == True:
# feature.qualifiers フィールドにD-loopの記載がある場合があるため
if 'note' in feature.qualifiers:
note = feature.qualifiers['note'][0]
else:
note = 'No Annotation'
sequence_slice = feature.location.extract(parent_sequence)
# 配列の末端
strand_int = feature.location.strand
if strand_int == 1:
strand = 'top' # 5'末端
elif strand_int == -1:
strand = 'btm' # 3'末端
else:
strand = 'cannotbe'
# 配列の開始位置
start = feature.location.parts[0].start.position + 1
# 配列の終了位置
end = feature.location.parts[-1].end.position
# 転写産物？の情報
product = 'NA'
# Fasta形式にした時の情報部分の作成
id_line = acc + '_' + organism + '_' + note + '_'
id_line += ' (' + str(start) + '..' + str(end) + ', strand=' + strand + ')'
# レコードから取得した情報をFasta形式にする
record_list = makeSeqRecord(sequence_slice, id_line)
# Fastaファイルの出力
fasta_file_name = './Results/Fasta/' + acc +'_' + feature.type + '.fasta'
type_count = len(record_list)
with open(fasta_file_name, mode = 'w') as ffn:
if type_count > 0:
SeqIO.write(record_list, ffn, 'fasta')
else :
ffn.write('No Annotation Found')
# merge fasta file name
mfa = './Results/Fasta/mergefa.fas'
with open(mfa, mode = 'a') as mfa:
if type_count > 0:
SeqIO.write(record_list, mfa, 'fasta')
else:
pass
# リストに追加された種名を記入
add_organism_list = []
add_organism_list = add_organism_list.append(organism)
print(add_organism_list)
with open("./Results/list.txt", mode = 'w') as spl:
spl.write(organism + '\n')
else:
pass
DNA領域に構わず、対象IDの配列を一括取得する
# --- import module --- #
import subprocess as sp
import pandas as pd
import pathlib, glob, re, pprint, os, datetime, shutil, time, sys
from Bio import Entrez, SeqIO
# from Bio.GenBank.Record import Record
from Bio.SeqFeature import SeqFeature, FeatureLocation
from Bio.SeqRecord import SeqRecord
from datetime import timedelta
from time import sleep
from dateutil.relativedelta import relativedelta
from subprocess import Popen, PIPE
from tqdm import tqdm
from termcolor import colored, cprint
# --- フォルダーのチェック --- #
def check_results_folder():
# Reusltsフォルダの有無のチェック
if pathlib.Path("./Results").is_dir() == True:
shutil.rmtree("./Results") # フォルダを消して
pathlib.Path("./Results").mkdir(exist_ok=True) # 作成
else:
pathlib.Path("./Results").mkdir(exist_ok=True)
# IDをまとめるファイルのチェック
if pathlib.Path('./Results/ID.txt').is_file == True:
pathlib.Path("./Results/ID.txt").unlink()
else:
with open('./Results/ID.txt', mode='w') as i:
i.write("")
# Fastaファイルをまとめるファイルのチェック
if pathlib.Path("./Results/Fasta").is_dir() == False:
pathlib.Path("./Results/Fasta").mkdir(exist_ok=True)
else:
pass
# --- チェック --- #
check_results_folder()
# --- Color setting --- #
class Color:
BLACK = '\033[30m'
RED = '\033[31m'
GREEN = '\033[32m'
YELLOW = '\033[33m'
BLUE = '\033[34m'
PURPLE = '\033[35m'
CYAN = '\033[36m'
WHITE = '\033[37m'
END = '\033[0m'
BOLD = '\038[1m'
UNDERLINE = '\033[4m'
INVISIBLE = '\033[08m'
REVERCE = '\033[07m'
# --- setting.txtの読み込み --- #
with open('Setting.txt', encoding = "utf-8") as st:
txt = st.read()
# --- NCBIへの自己紹介 --- #
Entrez.email = re.search(r'Email\s*=\s*(\S+)', txt).group(1)
Entrez.api_key = re.search(r'Api_key\s*=\s*(\S+)', txt).group(1)
# --- setting.txtから変数の情報を取得する --- #
# データベースを取得
db = re.search(r'Database1\s*=\s*(\S+)', txt).group(1)
# ファイルをダウンロードするときの形式
rettype = re.search(r'Rettype1\s*=\s*(\S+)', txt).group(1)
# データをダウンロードするときのフォーマット
retmode = re.search(r'Retmode1\s*=\s*(\S+)', txt).group(1)
# 1回の要求でダウンロードするデータ数
count_in_query = re.search(r'Count_in_query1\s*=\s*(\S+)', txt).group(1)
# いつからのデータを取得するか
minday = re.search(r'minday\s*=\s*(\S+)', txt).group(1)
# いつまでのデータを取得するか (常に今日になるようになっています)
maxday = datetime.date.today().strftime("%Y/%m/%d")
# 対象種を取得
species = re.search(r'Species\s*=\s*(.*)', txt).group(1)
# 対象遺伝子を設定
genetic= re.search(r'Genetic\s*=\s*(\S+)', txt).group(1)
# --- IDを検索する --- #
handle_search = Entrez.esearch(db = db, retmax = count_in_query, term = f'{species}[Organism] AND {genetic}[filter]')
record_search = Entrez.read(handle_search)
# 取得したIDのList
SeqID = record_search["IdList"]
# 操作しやすいようにID毎に改行する
str_ID = '\n'.join(SeqID)
# ResultsフォルダのID.txtに取得したIDを書き出す
with open('./Results/ID.txt', mode = 'a') as f:
f.write(str_ID + "\n")
# Dataframe の形でIDの読み込み
ID_dataframe = pd.read_csv("./Results/ID.txt", encoding="UTF-8", header=None)
# New function define
def makeSeqRecord(seq_obj, id_for_fasta):
record = SeqRecord(seq_obj, id=id_for_fasta, description="")
return record
# IDが記載されているID.txtの読み込み
ID_dataframe = pd.read_csv("./Results/ID.txt", encoding="UTF-8", header=None)
# Genbank ID数を取得
id_num = len(ID_dataframe)
# GWnbank ID の数だけForで情報を取得する
for i in tqdm(range(0, id_num), desc = Color.CYAN + 'Progress rate' + Color.END):
handle = Entrez.efetch(db = db, id = ID_dataframe.iloc[i,0], rettype = 'gb', retmode = 'text')
for record_info in SeqIO.parse(handle, "genbank"): # gbを指定するとGenbank形式の情報を取得する
# リスト格納
record_list = []
# 配列全体を格納
parent_sequence = record_info.seq
# 学名を格納
organism = record_info.annotations['organism']
# Accession numberを格納
acc = record_info.id
# 各領域ごとに配列を切り分けて配分
for feature in record_info.features:
# 不明な遺伝子領域 misc_feature 以外の場合
if feature.type == 'source':
sequence_slice = feature.location.extract(parent_sequence)
# 配列の末端取得
strand_int = feature.location.strand
if strand_int == 1:
strand = 'top' # 5'末端
elif strand_int == -1:
strand = 'btm' # 3'末端
else:
strand = 'cannotbe'
# 配列の開始位置指定
start = feature.location.parts[0].start.position + 1
# 配列の終了位置指定
end = feature.location.parts[-1].end.position
# Fasta形式にした時の情報部分(> ~~~ の部分)の作成
id_line = acc + '_' + organism + '(' + str(start) + '..' + str(end) + ', strand=' + strand + ')'
# レコードから取得した情報をFasta形式にする
record_list = makeSeqRecord(parent_sequence, id_line)
# アクセッションナンバー.fa のファスタファイルでファイル名を指定
fasta_file_name = './Results/Fasta/' + record_info.id + '.fasta'
# 一応0以上かチェック
type_count = len(record_list)
# 書き出し
with open(fasta_file_name, mode = 'w') as ffn:
if type_count > 0:
SeqIO.write(record_list, ffn, 'fasta')
else :
ffn.write('No Annotation Found')
# merge fasta file name
mfa = './Results/Fasta/mergefa.fas'
with open(mfa, mode = 'a') as mfa:
if type_count > 0:
SeqIO.write(record_list, mfa, 'fasta')
else:
pass
使用用途は何なのかと言われると難しい所ですが、誰かの仕事や研究の手助けとなれば幸いです。