Biopython入門 - 中編【NCBIデータベースやBlastへのアクセス】

2020年05月19日 2020年05月26日 12分26秒

Web master

この記事のタイトルとURLをコピーする

Biopythonはプログラミング言語Pythonをベースにした生命科学研究のためのパッケージで、その基礎的な使い方を3回に渡って紹介しています。

第1回目のBiopython入門 - 前編【インストール ~ 簡単な使い方】ではDNAやタンパク配列の作り方を中心に見てきました。

今回は、NCBIの代表的なツールであるEntrezデータベースとBlast検索をBiopythonで行う方法を紹介します。

この記事の内容

1 BiopythonからNCBI Entrezデータベースへのアクセス
2 BiopythonからNCBI Blastにアクセス
3 まとめに代えて
4 関連図書

BiopythonからNCBI Entrezデータベースへのアクセス

BiopythonからNCBIにアクセスする時には、まずEntrezユーザーガイドラインをよく読みましょう。

特に大事なところを2点抜粋して意訳すると、

Entrezユーザーガイドラインの特に重要な点

アクセスの頻度が一秒に三回以下になるようにする
自分のメールアドレスを指定する (問題があったときにNCBIから連絡ができるように)

この2点目のメールアドレスの指定については

from Bio import Entrez
Entrez.email = "[email protected]"

のように書けばOKです。

EInfoで利用可能なEntrezデータベースを取得する

NCBIのどのようなデータベースが利用可能なのかは、einfo ()を呼び出すことでXML形式で取得できます。

handle = Entrez.einfo()
result = handle.read()
print (result)

DbName (Database Name)で囲まれているデータベースが利用できるものになります。

EsearchでEntrezデータベースに検索する

利用可能なデータベースの一覧が取得できたので、そのうちの望みのものを選んで検索をかけてみます。

ここでは例としてPubmedにTP53についての文献がないか検索してみます (TP53はご存知の通り有名な腫瘍抑制遺伝子です)。

handle = Entrez.esearch(db="pubmed", term="TP53")
record = Entrez.read(handle)
record["IdList"]

するとこのように、TP53についての論文をPubmedで検索し、そのPMIDが返ってきます。

同様にして、ESearchを使ってGenBankデータベースに検索をかけてみます。ここでは同じくTP53遺伝子を対象として、人のデータに絞って(Homo Sapiens[Orgn]) やってみます。

handle = Entrez.esearch(db="nucleotide", term="Homo Sapiens[Orgn] AND TP53[Gene]")
record = Entrez.read(handle)
record["IdList"]

GenBankのID識別子は最後の["IdList"]で取得できました。

EFetchでデータをダウンロードする

取得したIDを使ったデータのダウンロードはEFetchを使って行うことができます。

例えば先ほど取得したTP53遺伝子, homo sapiensでの検索結果の一番最初にある1808862652について、そのデータを取得してみます。

handle = Entrez.efetch(db="nucleotide", id="1808862652", rettype="gb", retmode="text")
print(handle.read())

するとこのように、NCBI Nucleotideデータベースで検索したときと同じ形式でデータをダウンロードできました。

ここでは引数にrettype="gb"、retmode="text"を指定しているので、GenBankフォーマットでのダウンロードになりました (rettypeのgbはgenbankの略語です)。

もちろんFASTAフォーマットでもダウンロードでき、

rettype="fasta"

と変更するだけです。

handle = Entrez.efetch(db="nucleotide", id="1808862652", rettype="fasta", retmode="text")
print(handle.read())

それ以外のパラメーターの意味については、公式ドキュメントに紹介されています。

BiopythonからNCBI Blastにアクセス

BiopythonからNCBIの配列検索サービスBlastにアクセスするにはBio.Blastを使います。

Bio.Blastには、BLASTにアクセスするためのBio.Blast.NCBIWWWと、BLASTから返される検索結果 (XML形式) を見やすくしてくれるBio.Blast.NCBIXMLが含まれています。

例えば

from Bio.Blast import NCBIWWW
# 検索を起動する
# qblast (プログラム名、データベース名、検索キー)
result_handle = NCBIWWW.qblast ('blastn', 'nt', '8332116') 

# 検索結果を解釈する
from Bio.Blast import NCBIXML
blast_record = NCBIXML.read (result_handle)

BLASTのパラメーターは多数あり、qblastのパラメーターやBLASTのパラメーターの説明のページをご覧ください。

結果はblast_recordにあり、個々のアラインメントは.alignmentsで取得できます。

E_VALUE_THRESH = 0.05
for alignment in blast_record.alignments: # 個々のアラインメントを表示
  for hsp in alignment.hsps: 
    if hsp.expect < E_VALUE_THRESH: # もしE-valueが基準 (0.05) より低いなら
      print ('***Alignment***')
      print ('Sequence:', alignment.title) 
      print ('Length:', alignment.length) # アラインメント長
      print ('E value:', hsp.expect) # E-value
      print (hsp.query[0:75] + '.......') # 問い合わせ配列
      print (hsp.match[0:75] + '.......') # マッチしたか
      print (hsp.sbjct[0:75] + '.......') # ヒットした配列