BiopythonでNCBIWWWを用いてBlast検索を行う時

BiopythonのNCBIWWWを用いてBlast検索を実行することができます。
Blast検索自体は簡単だが、その後、XMLファイルの読み込みは煩雑であったため、備忘録として処理方法を記録する。

本ページではGoogle Colaboratoryを用いて実行したものになります。

https://biopython.org/docs/1.75/api/Bio.Blast.NCBIWWW.html

<h1>BLAST検索の実行</h1>
以下のコードにより、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"

データベースは以下のものに変更可能である。
<b>タンパク質データベースの場合</b>
"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)

<b>DNAデータベースの場合</b>
"nr" (nucleotide collection)
"wgs" (Whole-genome shotgun contigs)
"est" (Expressed sequence tag)
他

Format_typeは以下のものが選べる。
“HTML”, “Text”, “ASN.1”, “XML”。
特に指定しなければ、"XML"が選択される。

<b>query配列に関して</b>
input_seqを"seq.fasta"のようにfastaファイルの名前を指定して実行することも可能である。

<h2>その他追加で指定できるパラメーター</h2>

入力例
```
expect=10.0 #リスト化するe-valueの上限値
hitlist_size=50 #出力する遺伝子の数

```

<h2>中身の確認</h2>

以下の方法で中身を確認できる。また、ファイルとして保存も可能である。
```
#中身の印刷
print(result.getvalue())

#保存する方法
with open('myfile.txt', 'w', encoding='UTF-8') as f:
  f.write(result.getvalue())

```

<h1>データの読み込み</h1>

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)
```

中身は以下の方法で確認する。
確認方法は少しややこしい。

<h2>XMLファイルを読み込んだファイルないの階層構造について</h2>
おそらく中身は以下のような階層構造になっている。
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
    以下略
```

<h2>descriptionの解説</h2>
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

```

<h2>alignmentの情報の確認</h2>

以下の方法で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

```

<h2>MultipleAlignmentの情報の確認</h2>

以下の方法で確認できる。
しかし、通常の状態だと出力されないので、用いることはないだろう。
multipleAlignmentだけ、subjectの階層を持たない。
```
for i in blast_records:
  print(i.multiple_alignment)

#出力例
None

```

<h1>csvファイルと同じものを作成する (blastp出力結果を用いる場合)</h1>

上記情報を用いて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)
```