BiopythonのNCBIWWWを用いてBlast検索を実行することができます。
Blast検索自体は簡単だが、その後、XMLファイルの読み込みは煩雑であったため、備忘録として処理方法を記録する。
本ページではGoogle Colaboratoryを用いて実行したものになります。
BLAST検索の実行
以下のコードにより、BLAST検索の実行が可能です。from Bio.Blast import NCBIXML
from Bio.Blast import NCBIWWW
#ファイルの処理に必要なモジュールも同時に読み込んでおきます。
import io
import csv
#解析するための配列の入力 (アミノ酸配列の場合)
input_seq = "MDMH....." #@param {type:"string"}
#blastpの実行
result = NCBIWWW.qblast("blastp", "nr", input_seq, format_type='XML')
#tblastnの実行
result = NCBIWWW.qblast("tblastn", "nr", input_seq, format_type='XML')
"blastp", "nr", input_seqがそれぞれ、検索方法、データベースの種類、入力配列を表し、これらが最低限必要である。
データベースはそれぞれ変更して検索可能である。
検索方法は以下の5つを選択できる
"blastn"
"blastp"
"blastx"
"tblastn"
"tblastx"
データベースは以下のものに変更可能である。
タンパク質データベースの場合
"nr" (non-redundant protein sequences)
"refseq_select" (Refseq selected proteins)
"refseq_protein" (reference proteins)
"landmark" (Model organisms)
"swissprot" (UniProtKB/Swiss-Prot)
"pataa" (Patented protein sequences)
"pdb" (protein data bank)
"env_nr" (Metagenomic proteins)
"tsa_nr" (Transcriptome shotgun assembly proteins)
DNAデータベースの場合
"nr" (nucleotide collection)
"wgs" (Whole-genome shotgun contigs)
"est" (Expressed sequence tag)
他
Format_typeは以下のものが選べる。
“HTML”, “Text”, “ASN.1”, “XML”。
特に指定しなければ、"XML"が選択される。
query配列に関して
input_seqを"seq.fasta"のようにfastaファイルの名前を指定して実行することも可能である。
その他追加で指定できるパラメーター
入力例
expect=10.0 #リスト化するe-valueの上限値
hitlist_size=50 #出力する遺伝子の数
中身の確認
以下の方法で中身を確認できる。また、ファイルとして保存も可能である。
#中身の印刷
print(result.getvalue())
#保存する方法
with open('myfile.txt', 'w', encoding='UTF-8') as f:
f.write(result.getvalue())
データの読み込み
XML形式の出力結果をNCBIXMLを用いて読み出すことができる。
https://biopython.org/docs/1.75/api/Bio.Blast.NCBIXML.html
#直接読み込む場合
blast_records = NCBIXML.parse(result_2.getvalue())
#一度保存したファイルから読み込む場合
with open('myfile.txt', 'w', encoding='UTF-8') as f:
blast_records = NCBIXML.parse(f)
中身は以下の方法で確認する。
確認方法は少しややこしい。
XMLファイルを読み込んだファイルないの階層構造について
おそらく中身は以下のような階層構造になっている。 query配列は複数の配列を同時に検索しても対応できるよう、複数配列の出力が格納されていることが前提になっているようである。 しかし、pythonの検索では複数配列の同時検索はできない。 以下のページを参考に記述しています。 https://bi.biopapyrus.jp/python/biopython/ncbixml.html queryのなかに、三つの階層が存在します （description, alignment, MultipleAlignment)。 description, alignmentの中にはそれぞれ、subjectに対応する情報が格納されています。query1
description
subject1
title
score
e
num_alignments
subject2
title
score
e
num_alignments
subject3
以下略
alignment
subject1
title
length
hsps
score
bits
expect
num_alignments
identities
positives
gaps
strand
frame
query
query_start
match
sbjct
sbjct_start
subject2
以下略
MultipleAlignment
alignment
query2
description
以下略
descriptionの解説
blast_records内に複数のquery由来の検索結果が格納されていることが前提になっている。 そのため、forを用いてquery情報ごとに分割する。そのご、情報を獲得する。#description情報を獲得する
for i in blast_records:
#query1の階層にアクセス
#次のコードでquery1のdescription情報を呼び出す。
for j in i.descriptions:
print(j)
#出力例
gb|KOV25761.1| histidine ammonia-lyase [Streptomyces sp. XY431] 2643.0 0.0
ref|WP_043468418.1| histidine ammonia-lyase [Kitasatospora sp. MBT66] 2640.0 0.0
gb|KJY39357.1| histidine ammonia-lyase [Streptomyces sp. NRRL S-495] 2639.0 0.0
ref|WP_030393440.1| histidine ammonia-lyase [Kitasatospora purpeofusca] 2637.0 0.0
print(j)をpring(j.title)などに置き換えるとその情報のみ得ることができる。
#description情報の出力例
for i in blast_records:
for j in i.descriptions:
print(j.title)
print(j.score)
print(j.e)
print(j.num_alignments)
#出力例
gb|KOV25761.1| histidine ammonia-lyase [Streptomyces sp. XY431] #title
2643.0 #score
0.0 #e
1 #num_alignments
ref|WP_043468418.1| histidine ammonia-lyase [Kitasatospora sp. MBT66]
2640.0
0.0
1
gb|KJY39357.1| histidine ammonia-lyase [Streptomyces sp. NRRL S-495]
2639.0
0.0
1
alignmentの情報の確認
以下の方法でalignmentsの情報を呼び出せる。
for i in blast_records:
for j in i.alignments:
print(j)
#出力例
#subject1
gb|KOV25761.1| histidine ammonia-lyase [Streptomyces sp. XY431]
Length = 527
#subject2
ref|WP_043468418.1| histidine ammonia-lyase [Kitasatospora sp. MBT66]
Length = 525
#subject3
gb|KJY39357.1| histidine ammonia-lyase [Streptomyces sp. NRRL S-495]
Length = 527
同様にprint(j)をprint(j.title)などに置き換えると追加情報を取得できる。
for i in blast_records:
for j in i.alignments:
print(j.title)
print(j.length)
print(j.hsps)
出力例
#subject1
gb|KOV25761.1| histidine ammonia-lyase [Streptomyces sp. XY431] #title
527 #length
[<Bio.Blast.Record.HSP object at 0x7fc7bbb96fd0>] #hsps
#subject2
ref|WP_043468418.1| histidine ammonia-lyase [Kitasatospora sp. MBT66]
525
[<Bio.Blast.Record.HSP object at 0x7fc7bbc19910>]
#subject3
gb|KJY39357.1| histidine ammonia-lyase [Streptomyces sp. NRRL S-495]
527
[<Bio.Blast.Record.HSP object at 0x7fc7bbc19d90>]
上記の方法ではhspsの中身を見ることはできない。
なので、hspsの中身を見る場合は以下のコードを使うと良い。
for i in blast_records:
for j in i.alignments:
for k in j.hsps:
print(k)
#出力例
#subject1
Score 2643 (1022 bits), expectation 0.0e+00, alignment length 527
Query: 1 MDMHSAAAAHAPLVQVGKADVGAEDVLAVARGNARVEIGPDAIAE...PLA 527
MDMHSAAAAHAPLVQVGKADVGAEDVLAVARGNARVEIGPDAIAE...PLA
Sbjct: 1 MDMHSAAAAHAPLVQVGKADVGAEDVLAVARGNARVEIGPDAIAE...PLA 527
#subject2
Score 2640 (1021 bits), expectation 0.0e+00, alignment length 525
Query: 3 MHSAAAAHAPLVQVGKADVGAEDVLAVARGNARVEIGPDAIAEMA...PLA 527
MHSAAAAHAPLVQVGKADVGAEDVLAVARGNARVEIGPDAIAEMA...PLA
Sbjct: 1 MHSAAAAHAPLVQVGKADVGAEDVLAVARGNARVEIGPDAIAEMA...PLA 525
#subject3
Score 2639 (1021 bits), expectation 0.0e+00, alignment length 527
Query: 1 MDMHSAAAAHAPLVQVGKADVGAEDVLAVARGNARVEIGPDAIAE...PLA 527
MDMHSAAAAHAPLVQVGKADVGAEDVLAVARGNARVEIGPDAIAE...PLA
Sbjct: 1 MDMHSAAAAHAPLVQVGKADVGAEDVLAVARGNARVEIGPDAIAE...PLA 527
以下の方法で、hsps内の欲しい情報のみを取得できる。
for i in blast_records:
for j in i.alignments:
for k in j.hsps:
print(k.score)
print(k.bits)
print(k.expect)
print(k.num_alignments)
print(k.identities)
print(k.positives)
print(k.strand)
print(k.frame)
print(k.query)
print(k.query_start)
print(k.match)
print(k.sbjct)
print(k.sbjct_start)
count += 1
#出力例
#subject1
2643.0 #score
1022.69 #bits
0.0 #expect
None #num_alignments
523 #identities
526 #positives
(None, None) #strand
(0, 0) #frame
MDMHSA•••• #query
1 #query_start
MDMHSA•••• #match
MDMHSA•••• #sbjct
1 #sbjct_start
#subject2
2640.0
1021.54
0.0
None
524
525
(None, None)
(0, 0)
MHSAAA••••
3
MHSAA••••
MHSAA••••
1
#subject3
2639.0
1021.15
0.0
None
522
526
(None, None)
(0, 0)
MDMH••••
1
MDMH••••
VL+GSG••••
MDMHS••••
1
MultipleAlignmentの情報の確認
以下の方法で確認できる。
しかし、通常の状態だと出力されないので、用いることはないだろう。
multipleAlignmentだけ、subjectの階層を持たない。
for i in blast_records:
print(i.multiple_alignment)
#出力例
None
csvファイルと同じものを作成する (blastp出力結果を用いる場合)
上記情報を用いてhit tableをcsv形式で作成する。
query acc.ver, subject acc.ver, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score, % positives
の順にデータが格納されている。
一度, XMLファイルを保存した後に、以下のコードでHit table用のファイルを作成することができる。
#myfile.txtと保存したXMLデータを読み込む。
with open("myfile.txt", mode = 'r', encoding = 'utf-8') as fh:
#XML形式のファイルを読み込む。
blast_records = NCBIXML.parse(fh)
count = 1
#query配列の名前を指定する。
query_name = "query1"
list_1 = [] #descriptions由来の情報を格納するlist
list_2 = [] #hsps由来の情報を格納するlist
for i in blast_records:
#descriptions由来の情報を収集
for j in i.descriptions:
temp_list = []
temp_list.append(query_name)
title = j.title
title_2 = title.split("|")
temp_list.append(title_2[1])
list_1.append(temp_list)
print(list_1)
#hsps由来の情報を収集
for j in i.alignments:
for k in j.hsps:
temp_list = []
temp_list.append(k.identities / k.align_length * 100) #identities
temp_list.append(k.align_length) #alignment length
temp_list.append(k.align_length - k.identities) #mismatches
temp_list.append(k.gaps) #gap opens
temp_list.append(k.query_start) #q. start
temp_list.append(k.query_end) #q. end
temp_list.append(k.sbjct_start) #s. start
temp_list.append(k.sbjct_end) #s. end
temp_list.append(k.expect) #e value
temp_list.append(k.bits) #bit score
temp_list.append(k.positives / k.align_length * 100) #% positives
list_2.append(temp_list)
print(list_2)
#list_1とlist_2を結合したlistの作成
list_3 = [] #結合したlistの格納先
#実際の結合動作
for i in range(0,len(list_1)):
temp_list = []
for j in list_1[i]:
temp_list.append(j)
for j in list_2[i]:
temp_list.append(j)
list_3.append(temp_list)
print(list_3)
#ファイルの保存
with open("out.csv", mode = 'w', encoding = 'utf-8') as f:
writer = csv.writer(f)
writer.writerows(list_3)