Biopython を利用したNCBIのEntrez データベースへのアクセス #Python

多少意訳したり冗長なところは省いたりしています。

過去にも翻訳を試みた方がいるようですが、放置されているようなので、改めて訳します。

誤字・誤訳の指摘、補足など大歓迎です。

以下翻訳

EntrezはPubMed,GenBank, GEO等のNCBIのデータベースに対する、ユーザー向けに作られたータ取得システムです。ブラウザから直接アクセスして手動でクエリを行うこともできますが、BiopythonのBio.Entrezモジュールを介したプログラムによるアクセスも可能です。後者はPythonプログラムからPubMedにアクセスしたりGenBankのレコードにアクセスしたりすることを可能にします。

Bio.EntrezはEntrez Programming Utilities(a.k.a EUtils)を利用しています。これは8個のツールからなり、詳細はNCBIのページからご覧になれます。それぞれのツールはBio.Entrez内の8つのPythonの関数と対応しており、以下ではそれを説明していきます。このモジュールはクエリ用のURLが正しいことを保証し、クエリの頻度が(NCBIの要求を満たすよう)最大でも一秒間に三度以下になることを保証します。

Entrez Programming Utilitiesによって返されるアウトプットは大抵XMLフォーマットで書かれています。パースの方法はいくつかあります。

Bio.Entrezのパーサーを使用してPythonのオブジェクトに変換する。
Python標準ライブラリのDOMパーサを用いる
Python標準ライブラリのSAX(Simple API for XML)パーサを用いる。
XMLのraw textを読み、文字列解析をして扱う

2と3についてはPythonのドキュメンテーションを参照してください。ここでは1. Bio.Entrezのパーサーについて説明していきます。

NCBIはDTD(Document Type Definition)ファイルをXMLファイルに保持された情報を説明するのに用いています。NCBIによって使用されるDTDファイルのほとんどはBiopythonのディストリビューションに含まれています。Bio.EntrezはこのDTDファイルを用いて取得したXMLのパースを行います。

まれに、特定のXMLファイルと対応するDTDファイルが、Biopythonのディストリビューションから欠けている場合があります。具体的には、NCBIがDTDをアップデートした際に発生する可能性があります。
このような場合、Entrez.readが警告メッセージを表示し、欠けているDTDファイルのURLと名前を提示します。そしてパーサはXMLのパースを続行するため、DTDファイルをweb上から取得しようと試みます。
しかし、DTDファイルがローカルにあった方がパースははるかに早く進みます。そのためこのような場合は、警告メッセージ中のURLからDTDをダウンロードし、DTD用のディレクトリ(...site-packages/Bio/Entrez/DTD)中においてください。書き込み権限がない場合は~/.biopython/Bio/Entrez/DTDにおいてください。こちらのホーム以下のディレクトリの方が先に読み込まれるため、新しいバージョンのDTDファイルを使用してほしい場合にもここにおくことができます。
また、(訳注: pipやcondaではなく)ソースからBiopythonをインストールした場合、ソースコード中の..site-packages/Bio/Entrez/DTDディレクトリにDTDを配置してBioPythonを再インストールすることでアップデートすることも可能です。

EUtilsは出力を別のフォーマットにすることもできます。例えばFastaやシーケンスデータベース用のGenBank等です。あるいは文献データベースのMedLine用等です。詳しくは9.12説で見ていきます。

9.1 Entrezのガイドライン

Biopythonを用いてNCBIにアクセスする前に、NCBIのEntrezユーザーガイドラインを読んでください。NCBIがシステムの悪用を検出した場合、そのユーザーはアクセスを禁止される可能性があります！

要約すると

100以上のリクエストを投げる場合は、週末かUSAのピークタイム以外の時間に行いましょう。これは努力目標です。
NCBIの通常のアドレスではなく http://eutils.ncbi.nlm.nih.gov/ を用いましょう。Biopythonはこれを用いています。
アクセスの頻度が一秒に三度以下になるようにしましょう。(2009年以前は一秒に一度以下だったので、多少緩めになりました。)これはBiopythonが自動的に強制します。
オプションとして自身のemailアドレスを設定できます。過去には任意でしたが、ぜひ指定しましょう。これは問題が発生した時にNCBIがあなたに連絡できるようにするための設定です。呼び出しのたびに関数に引数として渡すか、以下のようにグローバルに設定するかの2通りのやり方があります。

from Bio import Entrez
Entrez.email = "A.N.Other@example.com"

こうしておけばBio.EntrezはEntrezへの呼び出しのたびにこのメールアドレスを用います。"A.N.Other@example.com"はドキュメンテーション用に特別に設けられたアドレスです(RFC 2606)。ランダムなメールアドレスを用いるのは絶対にやめてください。それよりは設定しない方がまだましです。2010年6月以降はこの設定が必須になりました。使用量が限度を超えた場合、E-utilitiesのアクセスをブロックする前に、NCBIはこのメールアドレスを介してユーザーにコンタクトを試みます。
* より大きなソフトウェアの一部としてBiopythonを使用したい場合、ツールパラメータを指定しましょう。これは関数呼び出しのたびに引数として与えるか(e.g. tool="MyLocalScript"を引数に加える。)、以下のようにグローバルにツール名を設定することができます。

from Bio import Entrez
Entrez.tool = "MyLocalScript"

ツール名のデフォルトはBiopythonです。

大規模なクエリを行う場合、NCBIはセッション履歴機能を用いることを推奨しています。(これにはWebEnvセッションクッキー文字列を使用します。詳しくは9.15節を見てください。)これによりほんの少しだけプログラムが複雑になります。

結論は、あなたの使用レベルに応じた対応を取ってください、ということになります。大量のデータのダウンロードを考えている場合、別の方法を検討してください。例えば、ヒトの全遺伝子への簡単なアクセス法が欲しい場合、それぞれの染色体上のものをGenBankファイルとしてFTPでダウンロードし、自身のBioSQLデータベース(20.5節参照)にインポートしましょう。

EInfo: Entrezデータベースの情報の取得

Einfoはフィールドのインデックスカウント(訳注: ?)、最終アップデート時間、それぞれのNCBIデータべースへアクセス可能なリンクを提供します。
加えて、Entrezユーティリティからアクセス可能な全データベースの名前を取得することもできます。


from Bio import Entrez
Entrez.email = "A.N.Other@example.com" # NCBIに自己紹介を怠らないこと
handle = Entrez.einfo()
result = handle.read()

result変数には利用可能なデータベースのリストがXML形式で入っています。

print(result)


<?xml version="1.0"?>
<!DOCTYPE eInfoResult PUBLIC "-//NLM//DTD eInfoResult, 11 May 2002//EN"
 "http://www.ncbi.nlm.nih.gov/entrez/query/DTD/eInfo_020511.dtd">
<eInfoResult>
<DbList>
        <DbName>pubmed</DbName>
        <DbName>protein</DbName>
        <DbName>nucleotide</DbName>
        <DbName>nuccore</DbName>
        <DbName>nucgss</DbName>
        <DbName>nucest</DbName>
        <DbName>structure</DbName>
        <DbName>genome</DbName>
        <DbName>books</DbName>
        <DbName>cancerchromosomes</DbName>
        <DbName>cdd</DbName>
        <DbName>gap</DbName>
        <DbName>domains</DbName>
        <DbName>gene</DbName>
        <DbName>genomeprj</DbName>
        <DbName>gensat</DbName>
        <DbName>geo</DbName>
        <DbName>gds</DbName>
        <DbName>homologene</DbName>
        <DbName>journals</DbName>
        <DbName>mesh</DbName>
        <DbName>ncbisearch</DbName>
        <DbName>nlmcatalog</DbName>
        <DbName>omia</DbName>
        <DbName>omim</DbName>
        <DbName>pmc</DbName>
        <DbName>popset</DbName>
        <DbName>probe</DbName>
        <DbName>proteinclusters</DbName>
        <DbName>pcassay</DbName>
        <DbName>pccompound</DbName>
        <DbName>pcsubstance</DbName>
        <DbName>snp</DbName>
        <DbName>taxonomy</DbName>
        <DbName>toolkit</DbName>
        <DbName>unigene</DbName>
        <DbName>unists</DbName>
</DbList>
</eInfoResult>

割とシンプルなXMLファイルなので、文字列検索でも情報を抜き取ることができますが、ここではBio.Entrezのパーサーを用いてXMLをPythonのオブジェクトに変換します。


from Bio import Entrez
handle = Entrez.einfo()
record = Entrez.read(handle)

ここでrecordはキーが一つだけの辞書です。

>>> record.keys()
... [u'DbList']

このキーに対応する値は上で紹介したXML中のデータベース名のリストです。

>>> record["DbList"]
['pubmed', 'protein', 'nucleotide', 'nuccore', 'nucgss', 'nucest',
 'structure', 'genome', 'books', 'cancerchromosomes', 'cdd', 'gap',
 'domains', 'gene', 'genomeprj', 'gensat', 'geo', 'gds', 'homologene',
 'journals', 'mesh', 'ncbisearch', 'nlmcatalog', 'omia', 'omim', 'pmc',
 'popset', 'probe', 'proteinclusters', 'pcassay', 'pccompound',
 'pcsubstance', 'snp', 'taxonomy', 'toolkit', 'unigene', 'unists']

それぞれのデータベースについて、EInfoを再度用いて、より詳細な情報を得ることができます。

>>> handle = Entrez.einfo(db="pubmed")
>>> record = Entrez.read(handle)
>>> record["DbInfo"]["Description"]
'PubMed bibliographic record'
>>> record["DbInfo"]["Count"]
'17989604'
>>> record["DbInfo"]["LastUpdate"]
'2008/05/24 06:45'

record["DbInfo"].keys()を用いて、このレコード中の他の情報を見てみることをお勧めします。ESearchを用いて検索可能なフィールドの情報が特に有用です。


>>> for field in record["DbInfo"]["FieldList"]:
...     print("%(Name)s, %(FullName)s, %(Description)s" % field)
ALL, All Fields, All terms from all searchable fields
UID, UID, Unique number assigned to publication
FILT, Filter, Limits the records
TITL, Title, Words in title of publication
WORD, Text Word, Free text associated with publication
MESH, MeSH Terms, Medical Subject Headings assigned to publication
MAJR, MeSH Major Topic, MeSH terms of major importance to publication
AUTH, Author, Author(s) of publication
JOUR, Journal, Journal abbreviation of publication
AFFL, Affiliation, Author's institutional affiliation and address
...

かなり長いリストですね。ということはしかし、PubMedのデータベースに関しては、特定の著者を探すのにJones[AUTH]としたり、Sanger[AFFL]でSanger Center所属の著者のみを探すといったことができるということを意味しています。とても扱いやすいでしょう。特定のデータベースに関して詳しくない場合は特に。

9.3 Esearch: Entrezデータベースのサーチ

これらのデータベースをサーチする場合、Bio.Entrez.esearch()を用います。例として、PubMedからBiopythonに関係する論文を検索してみましょう。

>>> from Bio import Entrez
>>> Entrez.email = "A.N.Other@example.com"     # Always tell NCBI who you are
>>> handle = Entrez.esearch(db="pubmed", term="biopython")
>>> record = Entrez.read(handle)
>>> record["IdList"]
['19304878', '18606172', '16403221', '16377612', '14871861', '14630660', '12230038']

7つのPubMed IDが出力されています。(19304878がBiopythonアプリケーションノートのPMIDです。)EFetchでそれぞれ取得できます。(9.6節を参照)

ESearchはGenBankの検索にも用いることができます。ここではCypripediodeae亜科のラン植物のmatK遺伝子を検索してみましょう。(Einfoを用いて、Entrezデータベースから検索可能なフィールドを知る方法については、9.2節を参照)

>>> handle = Entrez.esearch(db="nucleotide", term="Cypripedioideae[Orgn] AND matK[Gene]")
>>> record = Entrez.read(handle)
>>> record["Count"]
'25'
>>> record["IdList"]
['126789333', '37222967', '37222966', '37222965', ..., '61585492']

それぞれのID(126789333, 37222967, 37222966, …)はGenBankの識別子です。これらのGenBank recordをダウンロードしたい場合は9.6節節を参照してください。

Cypripedioideae[Orgn]のようにのように生物種名を用いる代わりにNCBIのtaxon identifierを用いてクエリ対象を絞り込むことができます。今回の場合これはtxid158330[Orgn]となります。今のところ、これはESearchのヘルプページにはドキュメントされていません。email付きのクエリに対してNCBIが回答してくれたことから判明しています。EntrezのWeb Interfaceで遊んでみると、検索文字列のフォーマットがわかります。例えばゲノムを検索する際にcomplete[prop]と入力すると、完全なゲノムのみが対象となります。

最後の例として、計算機科学関連のジャーナルのタイトルを取得しましょう。

>>> handle = Entrez.esearch(db="nlmcatalog", term="computational[Journal]", retmax='20')
>>> record = Entrez.read(handle)
>>> print("{} computational journals found".format(record["Count"]))
117 computational Journals found
>>> print("The first 20 are\n{}".format(record['IdList']))
['101660833', '101664671', '101661657', '101659814', '101657941',
 '101653734', '101669877', '101649614', '101647835', '101639023',
 '101627224', '101647801', '101589678', '101585369', '101645372',
 '101586429', '101582229', '101574747', '101564639', '101671907']

個々のジャーナルのIDからより詳細な情報を得たい際はやはりEFetchを用います。

Esearchにはほかにも有用なオプションが数多くあります。詳しくはヘルプページを見ましょう

EPost: IDのアップロード

EPostはその後の検索に使用するUIのリストをアップロードします。詳しくはEpostのヘルプページを見ましょう。
BiopythonではBio.Entrez.epost()から使用できます。

これが有用になる状況の例として、EFetchでダウンロードしたいIDの長いリストが手元にある状況を考えましょう。(塩基配列かもしれませんし、サイテーションかもしれません。なんでもよいです。)EFetchでそれらのリストの内容をリクエストした際、一つの長いURLが作られ、サーバーに送られます。リストがあまりにも長すぎると、URLが長すぎるためクエリが壊れる可能性があります。(例えばプロキシサーバがうまく対処できない可能性があります。)

対処法として、作業を2つのステップに分解するのが有用かもしれません。まずIDをEPostでアップロードし(内部ではHTMLのGETメソッドではなくPOSTメソッドを使用しているため、URLが長すぎる問題を回避できます。)向こうのサーバーが履歴(History)機能を持っているため、この長いIDのリストを参照してEFetchすることができます。

EPostがどのように機能するか、簡単な例とともに見ていきましょう。PuBMedのIDをアップロードします。


>>> from Bio import Entrez
>>> Entrez.email = "A.N.Other@example.com"     # Always tell NCBI who you are
>>> id_list = ["19304878", "18606172", "16403221", "16377612", "14871861", "14630660"]
>>> print(Entrez.epost("pubmed", id=",".join(id_list)).read())
<?xml version="1.0"?>
<!DOCTYPE ePostResult PUBLIC "-//NLM//DTD ePostResult, 11 May 2002//EN"
 "http://www.ncbi.nlm.nih.gov/entrez/query/DTD/ePost_020511.dtd">
<ePostResult>
 <QueryKey>1</QueryKey>
 <WebEnv>NCID_01_206841095_130.14.22.101_9001_1242061629</WebEnv>
</ePostResult>

返り値のXMLには二つの重要なフィールドが含まれています。QueryKeyとWebEnvです。両方を合わせることで、向こう側にあるあなたの履歴を参照できます。それぞれの値を別のEntrez呼び出しに使用しましょう。ここではEFetchです。

>>> from Bio import Entrez
>>> Entrez.email = "A.N.Other@example.com"     # Always tell NCBI who you are
>>> id_list = ["19304878", "18606172", "16403221", "16377612", "14871861", "14630660"]
>>> search_results = Entrez.read(Entrez.epost("pubmed", id=",".join(id_list)))
>>> webenv = search_results["WebEnv"]
>>> query_key = search_results["QueryKey"]

9.15節でこの履歴機能の詳しい使い方を見ていきます。

ESummary: primary IDからの要約の取得

ESummaryは文書の要約をIDから取得します。(より詳しくはHelpを参照してください。)上で検索した結果を用いて、例えばID30367のジャーナルを詳しい内容を見ていきましょう。

>>> from Bio import Entrez
>>> Entrez.email = "A.N.Other@example.com"     # Always tell NCBI who you are
>>> handle = Entrez.esummary(db="nlmcatalog", id="101660833")
>>> record = Entrez.read(handle)
>>> info = record[0]['TitleMainList'][0]
>>> print("Journal info\nid: {}\nTitle: {}".format(record[0]["Id"], info["Title"]))
Journal info
id: 101660833
Title: IEEE transactions on computational imaging.

EFetch: Entrezからの完全な情報のダウンロード

EFetchは全レコードをEntrezからダウンロードしたい際に利用するものです。ダウンロード可能なデータベースは複数あり、詳しくはHelpに書いてあります。

大抵のNCBIのデータベースは複数の異なるフォーマットをカバーしています。特定のファイルフォーマットを指定してダウンロードしたい場合はBio.Entrez.efetch()にrettypeあるいはretmodeという引数を与えてやる必要があります。それぞれのデータベースに応じた組み合わせはNCBIのefetchのウェブページにあります。(e.g. literature, sequences and taxonomy).

よく使われる例としては、

FASTAで配列をダウンロード
GenBank/GenPeptのプレインテキストフォーマットをダウンロード

等があります。(Bio.SeqIOでパースできます。詳しくは5.3.1節、9.6節を参照してください。)上で用いたCypripediodeae亜科の例から、GenBankのレコード番号186972394の情報をBio.Entrez.efetchを用いて取得しましょう。

>>> from Bio import Entrez
>>> Entrez.email = "A.N.Other@example.com"     # Always tell NCBI who you are
>>> handle = Entrez.efetch(db="nucleotide", id="186972394", rettype="gb", retmode="text")
>>> print(handle.read())
LOCUS       EU490707                1302 bp    DNA     linear   PLN 05-MAY-2008
DEFINITION  Selenipedium aequinoctiale maturase K (matK) gene, partial cds;
            chloroplast.
ACCESSION   EU490707
VERSION     EU490707.1  GI:186972394
KEYWORDS    .
SOURCE      chloroplast Selenipedium aequinoctiale
  ORGANISM  Selenipedium aequinoctiale
            Eukaryota; Viridiplantae; Streptophyta; Embryophyta; Tracheophyta;
            Spermatophyta; Magnoliophyta; Liliopsida; Asparagales; Orchidaceae;
            Cypripedioideae; Selenipedium.
REFERENCE   1  (bases 1 to 1302)
  AUTHORS   Neubig,K.M., Whitten,W.M., Carlsward,B.S., Blanco,M.A.,
            Endara,C.L., Williams,N.H. and Moore,M.J.
  TITLE     Phylogenetic utility of ycf1 in orchids
  JOURNAL   Unpublished
REFERENCE   2  (bases 1 to 1302)
  AUTHORS   Neubig,K.M., Whitten,W.M., Carlsward,B.S., Blanco,M.A.,
            Endara,C.L., Williams,N.H. and Moore,M.J.
  TITLE     Direct Submission
  JOURNAL   Submitted (14-FEB-2008) Department of Botany, University of
            Florida, 220 Bartram Hall, Gainesville, FL 32611-8526, USA
FEATURES             Location/Qualifiers
     source          1..1302
                     /organism="Selenipedium aequinoctiale"
                     /organelle="plastid:chloroplast"
                     /mol_type="genomic DNA"
                     /specimen_voucher="FLAS:Blanco 2475"
                     /db_xref="taxon:256374"
     gene            <1..>1302
                     /gene="matK"
     CDS             <1..>1302
                     /gene="matK"
                     /codon_start=1
                     /transl_table=11
                     /product="maturase K"
                     /protein_id="ACC99456.1"
                     /db_xref="GI:186972395"
                     /translation="IFYEPVEIFGYDNKSSLVLVKRLITRMYQQNFLISSVNDSNQKG
                     FWGHKHFFSSHFSSQMVSEGFGVILEIPFSSQLVSSLEEKKIPKYQNLRSIHSIFPFL
                     EDKFLHLNYVSDLLIPHPIHLEILVQILQCRIKDVPSLHLLRLLFHEYHNLNSLITSK
                     KFIYAFSKRKKRFLWLLYNSYVYECEYLFQFLRKQSSYLRSTSSGVFLERTHLYVKIE
                     HLLVVCCNSFQRILCFLKDPFMHYVRYQGKAILASKGTLILMKKWKFHLVNFWQSYFH
                     FWSQPYRIHIKQLSNYSFSFLGYFSSVLENHLVVRNQMLENSFIINLLTKKFDTIAPV
                     ISLIGSLSKAQFCTVLGHPISKPIWTDFSDSDILDRFCRICRNLCRYHSGSSKKQVLY
                     RIKYILRLSCARTLARKHKSTVRTFMRRLGSGLLEEFFMEEE"
ORIGIN
        1 attttttacg aacctgtgga aatttttggt tatgacaata aatctagttt agtacttgtg
       61 aaacgtttaa ttactcgaat gtatcaacag aattttttga tttcttcggt taatgattct
      121 aaccaaaaag gattttgggg gcacaagcat tttttttctt ctcatttttc ttctcaaatg
      181 gtatcagaag gttttggagt cattctggaa attccattct cgtcgcaatt agtatcttct
      241 cttgaagaaa aaaaaatacc aaaatatcag aatttacgat ctattcattc aatatttccc
      301 tttttagaag acaaattttt acatttgaat tatgtgtcag atctactaat accccatccc
      361 atccatctgg aaatcttggt tcaaatcctt caatgccgga tcaaggatgt tccttctttg
      421 catttattgc gattgctttt ccacgaatat cataatttga atagtctcat tacttcaaag
      481 aaattcattt acgccttttc aaaaagaaag aaaagattcc tttggttact atataattct
      541 tatgtatatg aatgcgaata tctattccag tttcttcgta aacagtcttc ttatttacga
      601 tcaacatctt ctggagtctt tcttgagcga acacatttat atgtaaaaat agaacatctt
      661 ctagtagtgt gttgtaattc ttttcagagg atcctatgct ttctcaagga tcctttcatg
      721 cattatgttc gatatcaagg aaaagcaatt ctggcttcaa agggaactct tattctgatg
      781 aagaaatgga aatttcatct tgtgaatttt tggcaatctt attttcactt ttggtctcaa
      841 ccgtatagga ttcatataaa gcaattatcc aactattcct tctcttttct ggggtatttt
      901 tcaagtgtac tagaaaatca tttggtagta agaaatcaaa tgctagagaa ttcatttata
      961 ataaatcttc tgactaagaa attcgatacc atagccccag ttatttctct tattggatca
     1021 ttgtcgaaag ctcaattttg tactgtattg ggtcatccta ttagtaaacc gatctggacc
     1081 gatttctcgg attctgatat tcttgatcga ttttgccgga tatgtagaaa tctttgtcgt
     1141 tatcacagcg gatcctcaaa aaaacaggtt ttgtatcgta taaaatatat acttcgactt
     1201 tcgtgtgcta gaactttggc acggaaacat aaaagtacag tacgcacttt tatgcgaaga
     1261 ttaggttcgg gattattaga agaattcttt atggaagaag aa
//

引数にrettype="gb"、retmode="text"を指定しているため、GenBankフォーマットでのダウンロードになっています。

2009年のイースターの日まで、EntrezのEFetch APIは"genbank"を返り値としていましたが、現在は"gb"、"gbwithparts"、(あるいはたんぱく質の場合"gp")を推奨するという旨の説明がオンラインにあります。また、2012年までデフォルトでプレーンテキストを返していましたが、現在はXMLを返しています。

上の例の代わりにrettype="fasta"でfastaフォーマットを取得することもできます。他のオプションはEFetch　Sequencesのヘルプページを見てください。
オプションの種類は、クエリ先のデータベースによるということを忘れないように！EFetchのヘルプを見ましょう。

Bio.SeqIOでパース可能なフォーマットでデータをダウンロードしたならば、SeqRecordですぐにパースできます。

>>> from Bio import Entrez, SeqIO
>>> handle = Entrez.efetch(db="nucleotide", id="186972394", rettype="gb", retmode="text")
>>> record = SeqIO.read(handle, "genbank")
>>> handle.close()
>>> print(record)
ID: EU490707.1
Name: EU490707
Description: Selenipedium aequinoctiale maturase K (matK) gene, partial cds; chloroplast.
Number of features: 3
...
Seq('ATTTTTTACGAACCTGTGGAAATTTTTGGTTATGACAATAAATCTAGTTTAGTA...GAA', IUPACAmbiguousDNA())

より典型的な利用法としては、配列をローカルファイルに保存し、その後Bio.seqIOにパースするといったケースがあげられるでしょう。こうした方が、何度もダウンロードしてNCBIのサーバーに負荷をかけるといったことがなくて済みます。

import os
from Bio import SeqIO
from Bio import Entrez
Entrez.email = "A.N.Other@example.com"  # Always tell NCBI who you are
filename = "gi_186972394.gbk"
if not os.path.isfile(filename):
    # Downloading...
    net_handle = Entrez.efetch(db="nucleotide", id="186972394", rettype="gb", retmode="text")
    out_handle = open(filename, "w")
    out_handle.write(net_handle.read())
    out_handle.close()
    net_handle.close()
    print("Saved")

print("Parsing...")
record = SeqIO.read(filename, "genbank")
print(record)

出力をXMLフォーマット(これならBio.Entrez.read()でパースできます。)で取得するにはretmode="xml"を使用しましょう。


>>> from Bio import Entrez
>>> handle = Entrez.efetch(db="nucleotide", id="186972394", retmode="xml")
>>> record = Entrez.read(handle)
>>> handle.close()
>>> record[0]["GBSeq_definition"]
'Selenipedium aequinoctiale maturase K (matK) gene, partial cds; chloroplast'
>>> record[0]["GBSeq_source"]
'chloroplast Selenipedium aequinoctiale'

Bio.Entrez.esearch()による検索を実行したのちに、Bio.Entrez.efetch()でダウンロードしたい場合、WebEnvの履歴機能を用いるべきです。9.15節を見ましょう。

ELink: NCBI Entrezの関連する要素を検索

ELinkをBiopythonから利用する場合はBio.Entrez.elink()を利用します。これはNCBI Entrezのデータベースから関連する要素を探すために用いられます。例えば遺伝子のデータベースから塩基配列のエントリを探すために使用できます。他にもcoolなことが色々できます。

ジャーナルBioinformaticsで2009年にパブリッシュされたBiopythonのアプリケーションノートに関連した論文を探すのにELinkを使用しましょう。このノートのPubMed IDは19304878ですので以下のようになります。

>>> from Bio import Entrez
>>> Entrez.email = "A.N.Other@example.com"
>>> pmid = "19304878"
>>> record = Entrez.read(Entrez.elink(dbfrom="pubmed", id=pmid))

recordはPythonのリストが入り、リストの要素が検索対象となったそれぞれのデータベースの情報に対応しています。今回は一つのPubMed IDのみを検索対象にしましたので、結果は1要素のリストです。この要素とは、我々の用いた検索タームおよびそれと関連する情報のうち見つかったものすべてを保持した辞書です。

>>> record[0]["DbFrom"]
'pubmed'
>>> record[0]["IdList"]
['19304878']

"LinkSetDB" keyの値には検索結果がリストとして入っています。リストのそれぞれの要素はそれぞれのデータベースに対応します。ここでは、全てPubMedデータベースにヒットしています。(もっとも、いくつかのサブカテゴリに分かれていますが。)

>>> len(record[0]["LinkSetDb"])
5
>>> for linksetdb in record[0]["LinkSetDb"]:
...     print(linksetdb["DbTo"], linksetdb["LinkName"], len(linksetdb["Link"]))
...
pubmed pubmed_pubmed 110
pubmed pubmed_pubmed_combined 6
pubmed pubmed_pubmed_five 6
pubmed pubmed_pubmed_reviews 5
pubmed pubmed_pubmed_reviews_five 5

実際の検索結果のデータは"Link" keyの値として保持されています。今回は110このヒットがあったようです。一つ目のヒットを見てみましょう。

>>> record[0]["LinkSetDb"][0]["Link"][0]
{u'Id': '19304878'}

これはクエリに用いたIDなので、ここでは特に役立ちません。2つ目の要素を見てみましょう。

>>> record[0]["LinkSetDb"][0]["Link"][1]
{u'Id': '14630660'}

このID14630660の論文はBiopythonのPDBパーサに関する物のようです。
forループで全てのPubMed IDをプリントしましょう。

>>> for link in record[0]["LinkSetDb"][0]["Link"]:
...     print(link["Id"])
19304878
14630660
18689808
17121776
16377612
12368254
......

うーん便利。しかし個人的には論文が引用されているかどうかが気になります。これは、少なくともPubmed Centralのジャーナルに関しては、ELinkで可能です。

ELinkそのものの詳細はELinkのヘルプを見ましょう。リンク名の説明だけのためにサブページが丸ごと一つ使われています。ここではデータベース間横断検索がどのように行われているかを説明しています。

EGQuery

EGQueryはGlobalなクエリを行い、検索ワードがそれぞれのEntrezデータベースで何回見つかるかを調べるます。それぞれのデータベースに大量のESearchをかけるよりも高速で便利です。以下の例ではBiopythonという用語のカウントを取得しています。

>>> from Bio import Entrez
>>> Entrez.email = "A.N.Other@example.com"     # Always tell NCBI who you are
>>> handle = Entrez.egquery(term="biopython")
>>> record = Entrez.read(handle)
>>> for row in record["eGQueryResult"]:
...     print(row["DbName"], row["Count"])
...
pubmed 6
pmc 62
journals 0
...

詳しくはEGQueryのヘルプを見てください。

ESpell: 特定の単語に関して推奨されたスペルを返す。

ここでは"Biopyhon"の推奨されるスペルを取得します。

>>> from Bio import Entrez
>>> Entrez.email = "A.N.Other@example.com"     # Always tell NCBI who you are
>>> handle = Entrez.espell(term="biopythooon")
>>> record = Entrez.read(handle)
>>> record["Query"]
'biopythooon'
>>> record["CorrectedQuery"]
'biopython'

より詳しくはヘルプを見てください。主な使用例としてはGUIツール作成時に自動的にsuggestionを出す、というのがあります。

Entrezの巨大なXMLファイルをパースする。

Entrez.read関数はEntrezから返されたXMLファイルをパースしてpythonオブジェクトに変換し、丸ごとメモリに載せます。メモリが足りなくなるような大きいファイルの場合はEntrez.parseを使用します。これはXMLを逐次的に読み込むジェネレータを返す関数です。

例えばEntrezの遺伝子データベースからftpで特定の生物種の遺伝子を丸ごとダウンロードしたとします。これはかなり巨大な可能性があります。実際2009年9月4日時点で、Entrezのヒト遺伝子データベースHomo_sapiens.ags.gzはASNフォーマットで116576kBのサイズでした。NCBIのgene2xmlというプログラムでXMLに変換することができます。(詳しくはNCBIのftpサイトを見てください)

gene2xml -b T -i Homo_sapiens.ags -o Homo_sapiens.xml

結果のXMLは6.1GBになります。Enrerz.readを実行するとたいていのパソコンではMemoryErrorがraiseされるでしょう。

このXMLファイルはそれぞれの遺伝子のEntrezのgeneレコードのリストです。Entrez.parseはこれを逐次的に読み込み、それぞれについてプリントしたり、情報を抜き出したりできます。例えば、以下はイテレートして遺伝子番号と名前をプリントします。

>>> from Bio import Entrez
>>> handle = open("Homo_sapiens.xml")
>>> records = Entrez.parse(handle)

>>> for record in records:
...     status = record['Entrezgene_track-info']['Gene-track']['Gene-track_status']
...     if status.attributes['value']=='discontinued':
...         continue
...     geneid = record['Entrezgene_track-info']['Gene-track']['Gene-track_geneid']
...     genename = record['Entrezgene_gene']['Gene-ref']['Gene-ref_locus']
...     print(geneid, genename)
This will print:

1 A1BG
2 A2M
3 A2MP
8 AA
9 NAT1
10 NAT2
11 AACP
12 SERPINA3
13 AADAC
14 AAMP
15 AANAT
16 AARS
17 AAVS1
...

9.11: エラーハンドリング

XMLのパース中に発生する問題として以下の3種類が考えられます。

そもそもXMLファイルでない
ファイルが途中で途切れている
XMLファイルの書式としては正しいが、セットとなるDTDに存在しない要素を持っている。

最初のケースの例としてfastaファイルをXMLとしてパースしようとしてしまった場合を見てみましょう

>>> from Bio import Entrez
>>> handle = open("NC_005816.fna") # a Fasta file
>>> record = Entrez.read(handle)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/local/lib/python2.7/site-packages/Bio/Entrez/__init__.py", line 257, in read
    record = handler.read(handle)
  File "/usr/local/lib/python2.7/site-packages/Bio/Entrez/Parser.py", line 164, in read
    raise NotXMLError(e)
Bio.Entrez.Parser.NotXMLError: Failed to parse the XML data (syntax error: line 1, column 0). Please make sure that the input data are in XM

ここでは、xmlファイルならば最初に必ず持っているはずの<?xml ...タグを見つけられなかったためにパーサが、これはXMLファイルではないと(正しく)判断しています。

2つ目の例の場合、パーサはCorrupetedXMLErrorをraiseします。以下は途中で途切れているxmlファイルの例です。

<?xml version="1.0"?>
<!DOCTYPE eInfoResult PUBLIC "-//NLM//DTD eInfoResult, 11 May 2002//EN" "http://www.ncbi.nlm.nih.gov/entrez/query/DTD/eInfo_020511.dtd">
<eInfoResult>
<DbList>
        <DbName>pubmed</DbName>
        <DbName>protein</DbName>
        <DbName>nucleotide</DbName>
        <DbName>nuccore</DbName>
        <DbName>nucgss</DbName>
        <DbName>nucest</DbName>
        <DbName>structure</DbName>
        <DbName>genome</DbName>
        <DbName>books</DbName>
        <DbName>cancerchromosomes</DbName>
        <DbName>cdd</DbName>

以下のようなトレースバックを吐きます。

>>> Entrez.read(handle)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/local/lib/python2.7/site-packages/Bio/Entrez/__init__.py", line 257, in read
    record = handler.read(handle)
  File "/usr/local/lib/python2.7/site-packages/Bio/Entrez/Parser.py", line 160, in read
    raise CorruptedXMLError(e)
Bio.Entrez.Parser.CorruptedXMLError: Failed to parse the XML data (no element found: line 16, column 0). Please make sure that the input data are not corrupted.

>>>

エラーメッセージを見れば、ファイルのどこでエラーが出たかがわかることに注意しましょう。

続いて3つ目の例です。例えば以下のxmlを見てください。

<?xml version="1.0"?>
<!DOCTYPE eInfoResult PUBLIC "-//NLM//DTD eInfoResult, 11 May 2002//EN" "http://www.ncbi.nlm.nih.gov/entrez/query/DTD/eInfo_020511.dtd">
<eInfoResult>
        <DbInfo>
        <DbName>pubmed</DbName>
        <MenuName>PubMed</MenuName>
        <Description>PubMed bibliographic record</Description>
        <Count>20161961</Count>
        <LastUpdate>2010/09/10 04:52</LastUpdate>
        <FieldList>
                <Field>
...
                </Field>
        </FieldList>
        <DocsumList>
                <Docsum>
                        <DsName>PubDate</DsName>
                        <DsType>4</DsType>
                        <DsTypeName>string</DsTypeName>
                </Docsum>
                <Docsum>
                        <DsName>EPubDate</DsName>
...
        </DbInfo>
</eInfoResult>

2行目でeInfo_020511.dtdがDTDファイルとして指定されていますが、何らかの理由で、<DocsumList>や他のいくつかのタグがこのDTDファイルに記載されていなかったとしましょう。その場合パーサーは途中でストップしValidationErrorをraiseします。

>>> from Bio import Entrez
>>> handle = open("einfo3.xml")
>>> record = Entrez.read(handle)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/local/lib/python2.7/site-packages/Bio/Entrez/__init__.py", line 257, in read
    record = handler.read(handle)
  File "/usr/local/lib/python2.7/site-packages/Bio/Entrez/Parser.py", line 154, in read
    self.parser.ParseFile(handle)
  File "/usr/local/lib/python2.7/site-packages/Bio/Entrez/Parser.py", line 246, in startElementHandler
    raise ValidationError(name)
Bio.Entrez.Parser.ValidationError: Failed to find tag 'DocsumList' in the DTD. To skip all tags that are not represented in the DTD, please call Bio.Entrez.read or Bio.Entrez.parse with validate=False.

そのようなタグを見つけた際、エラーを発生させずにスキップするようパーサーに指定することもできます。そのような場合、Entrez.read()あるいはEntrez.parse()にvalidate=Falseを引数として与えましょう。

>>> from Bio import Entrez
>>> handle = open("einfo3.xml")
>>> record = Entrez.read(handle, validate=False)
>>>

当然ながらDTDにはあってXMLにはない情報はパース結果に影響しない

9.12: 特別なパーサ

MedlineのパーサはBio.Medlineにあります。単一のMedlineのレコードのあるpubmed_result1.txtをパースしたいとしましょう。このファイルはBiopythonのTests\Medlineディレクトリに存在します。ファイルの中身は以下のようになっています。

PMID- 12230038
OWN - NLM
STAT- MEDLINE
DA  - 20020916
DCOM- 20030606
LR  - 20041117
PUBM- Print
IS  - 1467-5463 (Print)
VI  - 3
IP  - 3
DP  - 2002 Sep
TI  - The Bio* toolkits--a brief overview.
PG  - 296-302
AB  - Bioinformatics research is often difficult to do with commercial software. The
      Open Source BioPerl, BioPython and Biojava projects provide toolkits with
...

ファイルを開けてパースしましょう。

>>> from Bio import Medline
>>> with open("pubmed_result1.txt") as handle:
...    record = Medline.read(handle)
...

record変数にはMedlineのレコードがPythonの辞書として入っています。

>>> record["PMID"]
'12230038'
>>> record["AB"]
'Bioinformatics research is often difficult to do with commercial software.
The Open Source BioPerl, BioPython and Biojava projects provide toolkits with
multiple functionality that make it easier to create customised pipelines or
analysis. This review briefly compares the quirks of the underlying languages
and the functionality, documentation, utility and relative advantages of the
Bio counterparts, particularly from the point of view of the beginning
biologist programmer.'

ここでキーとして用いている値は少しあいまいです。

>>> help(record)

で簡単な要約を見ましょう。

複数のMedlineのレコードをパースしたい場合は、parse関数を使用しましょう。

>>> from Bio import Medline
>>> with open("pubmed_result2.txt") as handle:
...     for record in Medline.parse(handle):
...         print(record["TI"])
...
A high level interface to SCOP and ASTRAL implemented in python.
GenomeDiagram: a python package for the visualization of large-scale genomic data.
Open source clustering software.
PDB file parser and structure class implemented in Python.

ファイル中のMedlineレコードをパースする代わりにBio.Entrez.efetchで取得したMedlineのレコードをパースするkと思できます。PubMed中のBiopythonに関連した全MedLineレコードを見てみましょう。

>>> from Bio import Entrez
>>> Entrez.email = "A.N.Other@example.com"     # Always tell NCBI who you are
>>> handle = Entrez.esearch(db="pubmed", term="biopython")
>>> record = Entrez.read(handle)
>>> record["IdList"]
['19304878', '18606172', '16403221', '16377612', '14871861', '14630660', '12230038']

次いでBio.Entrez.efetchでMedlineをダウンロードします。

>>> idlist = record["IdList"]
>>> handle = Entrez.efetch(db="pubmed", id=idlist, rettype="medline", retmode="text")

ここではrettype="medline"を指定しています。ではこのレコードをBio.Medlineでパースしましょう。

>>> from Bio import Medline
>>> records = Medline.parse(handle)
>>> for record in records:
...     print(record["AU"])
['Cock PJ', 'Antao T', 'Chang JT', 'Chapman BA', 'Cox CJ', 'Dalke A', ..., 'de Hoon MJ']
['Munteanu CR', 'Gonzalez-Diaz H', 'Magalhaes AL']
['Casbon JA', 'Crooks GE', 'Saqi MA']
['Pritchard L', 'White JA', 'Birch PR', 'Toth IK']
['de Hoon MJ', 'Imoto S', 'Nolan J', 'Miyano S']
['Hamelryck T', 'Manderick B']
['Mangalam H']

比較のため、XMLフォーマットを用いた例を見ましょう。

>>> idlist = record["IdList"]
>>> handle = Entrez.efetch(db="pubmed", id=idlist, rettype="medline", retmode="xml")
>>> records = Entrez.read(handle)
>>> for record in records:
...     print(record["MedlineCitation"]["Article"]["ArticleTitle"])
Biopython: freely available Python tools for computational molecular biology and
 bioinformatics.
Enzymes/non-enzymes classification model complexity based on composition, sequence,
 3D and topological indices.
A high level interface to SCOP and ASTRAL implemented in python.
GenomeDiagram: a python package for the visualization of large-scale genomic data.
Open source clustering software.
PDB file parser and structure class implemented in Python.
The Bio* toolkits--a brief overview.

ここでは単純化のため、両例ともESearchとEFetchを組み合わせて使用していることに注目してください。
このような状況では、NCBI側は、9.15節で紹介している履歴機能を用いることを期待しているでしょう。

9.12.2

GEO(Gene Expression Omnibus)はNGSとマイクロアレイの発現量データのリポジトリです。GEOフォーマットのデータのパースにはBio.Geoモジュールが使用できます。

>>> from Bio import Geo
>>> handle = open("GSE16.txt")
>>> records = Geo.parse(handle)
>>> for record in records:
...     print(record)

"gds"データベースへの検索はESearchで行うことができます。

>>> from Bio import Entrez
>>> Entrez.email = "A.N.Other@example.com" # Always tell NCBI who you are
>>> handle = Entrez.esearch(db="gds", term="GSE16")
>>> record = Entrez.read(handle)
>>> record["Count"]
2
>>> record["IdList"]
['200000016', '100000028']

EntrezのWebサイトでは、UID "200000016"はGDS16に対応し、もう一つの"100000028"は対応するプラットフォームGPL28に対応します。残念ながらこの文書の執筆時点ではNCBIはEntrezを介したGEOファイルのダウンロードには対応していない用ですが、FTPでのダウンロードはこちらから簡単にできます ftp://ftp.ncbi.nih.gov/pub/geo/

訳注: 現在(2016/07/13)は、NCBI側は対応しているようですが、Biopythonはまだ対応していない用です。

今回の例で用いているデータのgzipされたものは ftp://ftp.ncbi.nih.gov/pub/geo/DATA/SOFT/by_series/GSE16/GSE16_family.soft.gz こちらからダウンロードできます。pythonでの解凍にはgzipモジュールを使用しましょう。

9.12.3 UniGeneのレコードのパース

UniGeneはNCBIのトランスクリプトームデータベースです。それぞれのUniGeneのレコードは特定の生物種の特定の遺伝子に対応する全転写産物です。典型的なUniGeneのレコードは以下のような見た目になっています。

ID          Hs.2
TITLE       N-acetyltransferase 2 (arylamine N-acetyltransferase)
GENE        NAT2
CYTOBAND    8p22
GENE_ID     10
LOCUSLINK   10
HOMOL       YES
EXPRESS      bone| connective tissue| intestine| liver| liver tumor| normal| soft tissue/muscle tissue tumor| adult
RESTR_EXPR   adult
CHROMOSOME  8
STS         ACC=PMC310725P3 UNISTS=272646
STS         ACC=WIAF-2120 UNISTS=44576
STS         ACC=G59899 UNISTS=137181
...
STS         ACC=GDB:187676 UNISTS=155563
PROTSIM     ORG=10090; PROTGI=6754794; PROTID=NP_035004.1; PCT=76.55; ALN=288
PROTSIM     ORG=9796; PROTGI=149742490; PROTID=XP_001487907.1; PCT=79.66; ALN=288
PROTSIM     ORG=9986; PROTGI=126722851; PROTID=NP_001075655.1; PCT=76.90; ALN=288
...
PROTSIM     ORG=9598; PROTGI=114619004; PROTID=XP_519631.2; PCT=98.28; ALN=288

SCOUNT      38
SEQUENCE    ACC=BC067218.1; NID=g45501306; PID=g45501307; SEQTYPE=mRNA
SEQUENCE    ACC=NM_000015.2; NID=g116295259; PID=g116295260; SEQTYPE=mRNA
SEQUENCE    ACC=D90042.1; NID=g219415; PID=g219416; SEQTYPE=mRNA
SEQUENCE    ACC=D90040.1; NID=g219411; PID=g219412; SEQTYPE=mRNA
SEQUENCE    ACC=BC015878.1; NID=g16198419; PID=g16198420; SEQTYPE=mRNA
SEQUENCE    ACC=CR407631.1; NID=g47115198; PID=g47115199; SEQTYPE=mRNA
SEQUENCE    ACC=BG569293.1; NID=g13576946; CLONE=IMAGE:4722596; END=5'; LID=6989; SEQTYPE=EST; TRACE=44157214
...
SEQUENCE    ACC=AU099534.1; NID=g13550663; CLONE=HSI08034; END=5'; LID=8800; SEQTYPE=EST
//

このレコードは、N-アセチルトランスフェラーゼをコードするヒトNAT2遺伝子由来の転写産物(SEQUENCE行)のセットを表しています。PROTSIM行はNAT2に有意に似ているたんぱく質を表しています。また、STS行はゲノム中の対応するsequence-tagged siteの場所を指しています。

UniGeneファイルのパースにはBio.UniGeneモジュールを使用しましょう。

>>> from Bio import UniGene
>>> input = open("myunigenefile.data")
>>> record = UniGene.read(input)

UniGene.readの返すレコードはPythonのオブジェクトでそのフィールドにはUniGeneのレコードのフィールドが入っています。例えば

>>> record.ID
"Hs.2"
>>> record.title
"N-acetyltransferase 2 (arylamine N-acetyltransferase)"

といった具合です。EXPRESSとRESTR_EXPR行はpythonのリストで、要素の型は文字列です。

['bone', 'connective tissue', 'intestine', 'liver', 'liver tumor', 'normal', 'soft tissue/muscle tissue tumor', 'adult']

STS、PROTSIM、SEQUENCE行には特別なオブジェクトが返されます。それぞれの行にある要素をフィールドとして持っています。

>>> record.sts[0].acc
'PMC310725P3'
>>> record.sts[0].unists
'272646'

PROTSIMやSEQUENCE行についても同様です。

一つ以上のUniGeneレコードを持つファイルをパースするにはBio.Unigeneのparse関数を使用しましょう。

>>> from Bio import UniGene
>>> input = open("unigenerecords.data")
>>> records = UniGene.parse(input)
>>> for record in records:
...     print(record.ID)

9.13: プロキシの使用

通常、プロキシに関して心を煩わせる必要はありません。しかし、あなたのネットワーク環境でプロキシサーバに関する問題が発生した場合、以下のように対処する必要があります。Bio.Entrezはurllibというpythonの標準ライブラリを用いてNCBIのサーバにアクセスします。これはhttp_proxyという環境変数を用いて、単純なプロキシならば自動でうまいこと対処してくれます。しかし残念なことに、このモジュールは認証の必要なプロキシには対処してくれません。

そのような場合の対処法は2つ

http_proxy環境変数をもう一度(訳注: 認証情報付きで)セットする。
Pythonスクリプトの最初で以下のように変数をセットする。

import os
os.environ["http_proxy"] = "http://proxyhost.example.com:8080"

詳しくはurllibのドキュメントを見てください。

9.14: 使用例

あなたが医療関係従事者であるか、ヒトを対象とした研究を行っているのであれば(あるいはそうでなくとも！)PubMedは情報源として実に有用です。他の場合と同様、pythonスクリプトで情報を抜き出してみましょう。

この例ではラン科植物と関係がある全ての論文を引き出すクエリをPubMedに投げてみます。(なぜランなのかは2.3節を参照)まずは論文の数を調べましょう。

>>> from Bio import Entrez
>>> Entrez.email = "A.N.Other@example.com"     # Always tell NCBI who you are
>>> handle = Entrez.egquery(term="orchid")
>>> record = Entrez.read(handle)
>>> for row in record["eGQueryResult"]:
...     if row["DbName"]=="pubmed":
...         print(row["Count"])
463

次にbio.Entrez.efetch関数でこれら463論文のPubmedIDを調べましょう。

>>> handle = Entrez.esearch(db="pubmed", term="orchid", retmax=463)
>>> record = Entrez.read(handle)
>>> idlist = record["IdList"]
>>> print(idlis)

これはPythonのPubMed IDの文字列を持つリストを返します。

['18680603', '18665331', '18661158', '18627489', '18627452', '18612381',
'18594007', '18591784', '18589523', '18579475', '18575811', '18575690',
...

取得したので、次は対応するMedlineのレコードと、それの持つ情報を引き出していきます。MedlineのフラットファイルフォーマットでダウンロードしBio.Medlineでパースします。

>>> from Bio import Medline
>>> handle = Entrez.efetch(db="pubmed", id=idlist, rettype="medline",
                           retmode="text")
>>> records = Medline.parse(handle)

ここでは検索とデータの取得を別々に行いましたが、NCBIは履歴機能を用いることを推奨しています。9.15節を参照。

レコードはイテレータであることを忘れずに。ですのでレコードの走査は一度だけ行えます。保持したい場合はリストに変換しましょう。

>>> records = list(records)

では、レコードをイテレートして、それぞれの情報をプリントしましょう。

>>> for record in records:
...     print("title:", record.get("TI", "?"))
...     print("authors:", record.get("AU", "?"))
...     print("source:", record.get("SO", "?"))
...     print("")

出力は以下のようになります。

title: Sex pheromone mimicry in the early spider orchid (ophrys sphegodes):
patterns of hydrocarbons as the key mechanism for pollination by sexual
deception [In Process Citation]
authors: ['Schiestl FP', 'Ayasse M', 'Paulus HF', 'Lofstedt C', 'Hansson BS',
'Ibarra F', 'Francke W']
source: J Comp Physiol [A] 2000 Jun;186(6):567-74

特に注目すべきはpythonのリストとして返された著者のリストです。これによりpythonの通常のツールを用いた検索が容易になります。例えば、特定の著者の大量のエントリを以下のようにしてループすることができます。

>>> search_author = "Waits T"

>>> for record in records:
...     if not "AU" in record:
...         continue
...     if search_author in record["AU"]:
...         print("Author %s found: %s" % (search_author, record["SO"]))

このセクションを読んで、MedlineとEntrezのインターフェイスの力と、どのように使用したらいいかアイディアが浮かべば幸いです。

9.14.2: Entrez Nucleotide レコードから検索、ダウンロード、パースする。

ここではリモートにあるEntrezに対してクエリを投げるシンプルな例を提示します。2.3節のパースの例では、NCBIのEntrezのウェブサイトを用いてNCBIの塩基配列データベースからユリ科植物(lady slipper orchids)の情報を取得する例を見てきました。ここではpythonスクリプトを用いてその操作を自動化する方法を見ていきます。この例ではEntrezモジュールのみで、接続し、結果を取得し、パースする方法を見ていきます。

まず、ダウンロードの前にEGQueryを用いてクエリがヒットする数のみを取得します。
今回、興味があるのは塩基配列のみです。

>>> from Bio import Entrez
>>> Entrez.email = "A.N.Other@example.com"     # Always tell NCBI who you are
>>> handle = Entrez.egquery(term="Cypripedioideae")
>>> record = Entrez.read(handle)
>>> for row in record["eGQueryResult"]:
...     if row["DbName"]=="nuccore":
...         print(row["Count"])
814

なるほど、つまり814の塩基配列がEntrezから手に入ることが期待されるわけですね。(これは2008年に公開されたバージョン1の情報で将来的にはおそらく増加すると思われます。)もしもヒット数がとんでもなく多い数になったとしたら、本当にすべてをダウンロードしたいのか再考した方が良いかもしれません。次のステップでそうします。

>>> from Bio import Entrez
>>> handle = Entrez.esearch(db="nucleotide", term="Cypripedioideae", retmax=814)
>>> record = Entrez.read(handle)

ここで、recordはpythonの辞書で、検索の結果がいくつかの追加情報とともに入っています。内容を見ていきましょう。

>>> print(record.keys())
[u'Count', u'RetMax', u'IdList', u'TranslationSet', u'RetStart', u'QueryTranslation']

まず、結果の個数をカウントしましょう。

>>> print(record["Count"])
'814'

期待した通りの個数です。

>>> len(record["IdList"])
814

では最初5つの結果を見てみましょう。

>>> record["IdList"][:5]
['187237168', '187372713', '187372690', '187372688', '187372686']

efetchを用いてこれらのレコードをダウンロードしましょう。一つ一つダウンロードすることも可能ですが、NCBI側の負担を減らすために、クエリは以下のように一度にまとめて投げた方が良いです。さらに言えば、9.15節で説明する履歴機能を用いるのが理想です。

>>> idlist = ",".join(record["IdList"][:5])
>>> print(idlist)
187237168,187372713,187372690,187372688,187372686
>>> handle = Entrez.efetch(db="nucleotide", id=idlist, retmode="xml")
>>> records = Entrez.read(handle)
>>> len(records)
5

これらそれぞれのレコードはGenBankのレコードに対応します。

>>> print(records[0].keys())
[u'GBSeq_moltype', u'GBSeq_source', u'GBSeq_sequence',
 u'GBSeq_primary-accession', u'GBSeq_definition', u'GBSeq_accession-version',
 u'GBSeq_topology', u'GBSeq_length', u'GBSeq_feature-table',
 u'GBSeq_create-date', u'GBSeq_other-seqids', u'GBSeq_division',
 u'GBSeq_taxonomy', u'GBSeq_references', u'GBSeq_update-date',
 u'GBSeq_organism', u'GBSeq_locus', u'GBSeq_strandedness']

>>> print(records[0]["GBSeq_primary-accession"])
DQ110336

>>> print(records[0]["GBSeq_other-seqids"])
['gb|DQ110336.1|', 'gi|187237168']

>>> print(records[0]["GBSeq_definition"])
Cypripedium calceolus voucher Davis 03-03 A maturase (matR) gene, partial cds;
mitochondrial

>>> print(records[0]["GBSeq_organism"])
Cypripedium calceolus

サッと検索したい場合はここで説明したようにするのが良いです。さらに重い処理をしたい場合、9.15節を見ましょう。

9.14.3 GenBankレコードの検索、ダウンロード、パース

GenBankレコードフォーマットは配列、配列の特徴、あるいは他の関連データの保存形式としてとてもポピュラーです。このフォーマットはNCBIのデータベース http://www.ncbi.nlm.nih.gov/ から情報を取得するのに推奨される方法です。

この例では、NCBIのデータベースにクエリを投げてレコードを取得し、Bio.SeqIO(5.3.1項で説明)を用いてパースする方法を見ていきます。単純化のためにWebEnvの履歴機能(9.15節で説明)は用いていません。

まず、クエリを行って取得したレコードのIDを探し出します。ここでは私が興味のあるサボテンの一種であるOpuntia(Pricly Pear Cactus)のクイックサーチを行い、対応する全GI(GenBank identifier)を取得します。まずはレコードの数を確認しましょう。


>>> from Bio import Entrez
>>> Entrez.email = "A.N.Other@example.com"     # Always tell NCBI who you are
>>> handle = Entrez.egquery(term="Opuntia AND rpl16")
>>> record = Entrez.read(handle)
>>> for row in record["eGQueryResult"]:
...     if row["DbName"]=="nuccore":
...         print(row["Count"])
...
9

次にGenBank IDのリストを取得します。

>>> handle = Entrez.esearch(db="nuccore", term="Opuntia AND rpl16")
>>> record = Entrez.read(handle)
>>> gi_list = record["IdList"]
>>> gi_list
['57240072', '57240071', '6273287', '6273291', '6273290', '6273289', '6273286',
'6273285', '6273284']

次にこれらのGIを用いてGenBankレコードを取得します。古いバージョンのBiopythonでは、GIの数字をカンマ区切りのリストにしてEntrezに渡す必要がありましたが、Biopython1.59以降ではリストを渡せば自動で変換してくれます。

>>> gi_str = ",".join(gi_list)
>>> handle = Entrez.efetch(db="nuccore", id=gi_str, rettype="gb", retmode="text")

もし生のGenBankファイルが見たければ、このhandleから読み出し、結果をプリントできます。

>>> text = handle.read()
>>> print(text)
LOCUS       AY851612                 892 bp    DNA     linear   PLN 10-APR-2007
DEFINITION  Opuntia subulata rpl16 gene, intron; chloroplast.
ACCESSION   AY851612
VERSION     AY851612.1  GI:57240072
KEYWORDS    .
SOURCE      chloroplast Austrocylindropuntia subulata
  ORGANISM  Austrocylindropuntia subulata
            Eukaryota; Viridiplantae; Streptophyta; Embryophyta; Tracheophyta;
            Spermatophyta; Magnoliophyta; eudicotyledons; core eudicotyledons;
            Caryophyllales; Cactaceae; Opuntioideae; Austrocylindropuntia.
REFERENCE   1  (bases 1 to 892)
  AUTHORS   Butterworth,C.A. and Wallace,R.S.
...

ここでは単に生のレコードを取得しています。もっとPythonに適した形で取得したければBio.SeqIOを用いてGenBankのデータをSeqFeatureのようなSeqRecordオブジェクトにパースすることができます。(詳しくは5章)

>>> from Bio import SeqIO
>>> handle = Entrez.efetch(db="nuccore", id=gi_str, rettype="gb", retmode="text")
>>> records = SeqIO.parse(handle, "gb")

レコードを走査して、興味のある情報を見ていくことができます。

>>> for record in records:
>>> ...    print("%s, length %i, with %i features" \
>>> ...           % (record.name, len(record), len(record.features)))
AY851612, length 892, with 3 features
AY851611, length 881, with 3 features
AF191661, length 895, with 3 features
AF191665, length 902, with 3 features
AF191664, length 899, with 3 features
AF191663, length 899, with 3 features
AF191660, length 893, with 3 features
AF191659, length 894, with 3 features
AF191658, length 896, with 3 features

このようにクエリを自動で行うことは手動で行うことに比べて大きな前進です。このモジュールはNCBIの「一秒にクエリは最大３回まで」というルールに従うとはいえ、NCBIが推奨しているルールには「ピークとなる時間を避けてクエリを行う」というようなものもあります(詳しくは9.1節)。単純化のため、ここではWebEnvの履歴機能(くわしくは9.15節)を使用していないことに注意してください。

最後に、もし解析を繰り返し行うつもりでいるのであれば、例にあるようにNCBIからダウンロードしたファイルを即座にパースするのではなくレコードをただ一度だけダウンロードしてハードディスクに保存し、ローカルファイルをパースするべきです。

9.14.4 生物種の系統樹を探す。

植物の例を用いた説明を続けていきます。アツモリソウ亜科(Cypripedioideae)の系統樹データベースを検索します。これはただ一つの系統樹IDを返します。

>>> from Bio import Entrez
>>> Entrez.email = "A.N.Other@example.com"     # Always tell NCBI who you are
>>> handle = Entrez.esearch(db="Taxonomy", term="Cypripedioideae")
>>> record = Entrez.read(handle)
>>> record["IdList"]
['158330']
>>> record["IdList"][0]
'158330'

では、efetchを用いて系統樹データベースのこのエントリをダウンロードし、パースします。

>>> handle = Entrez.efetch(db="Taxonomy", id="158330", retmode="xml")
>>> records = Entrez.read(handle)

このレコードもやはり、大量の情報を保持しています。

>>> records[0].keys()
[u'Lineage', u'Division', u'ParentTaxId', u'PubDate', u'LineageEx',
 u'CreateDate', u'TaxId', u'Rank', u'GeneticCode', u'ScientificName',
 u'MitoGeneticCode', u'UpdateDate']

このレコードから系統樹の情報を直接取得できます。

>>> records[0]["Lineage"]
'cellular organisms; Eukaryota; Viridiplantae; Streptophyta; Streptophytina;
 Embryophyta; Tracheophyta; Euphyllophyta; Spermatophyta; Magnoliophyta;
 Liliopsida; Asparagales; Orchidaceae'

このレコードデータには、ここで見たデータよりも更に大量の情報があります。例えばLineageではなく"LineageEx"以下のデータを見てください。NCBIの系統樹のIDも取得することができます。

9.15: 履歴機能(history)とWebEnv

9.15.1: 配列検索の検索とダウンロードのための履歴機能の使用

Opuntia(サボテンの一種)のrpl16の塩基配列を取得し、FASTAとして保存したいとしましょう。9.14.3節で紹介したように、単純にBio.Entrez.esearch()とBio.Entrez.efetch()を組み合わせて処理しすることができます。前者でリストを取得し、後者で塩基配列を全てダウンロードするといった具合です。

しかし、推奨されたアプローチはまず履歴機能を用いたサーチを行い、次いでリファレンスを用いて結果をダウンロードすることです。そのようにすればNCBI側でキャッシングなどの適切な対処を行うことができます。

このために、Bio.Entrez.esearch()を通常通り行います。ただし、余計な引数を一つ加えます。usehistory="y"です。

>>> from Bio import Entrez
>>> Entrez.email = "history.user@example.com"
>>> search_handle = Entrez.esearch(db="nucleotide",term="Opuntia[orgn] and rpl16",
                                   usehistory="y")
>>> search_results = Entrez.read(search_handle)
>>> search_handle.close()

こうしてXMLを取得しても、その中の検索結果に変わりはありません。

>>> gi_list = search_results["IdList"]
>>> count = int(search_results["Count"])
>>> assert count == len(gi_list)

しかし、２つの追加情報が得られます。WebEnvというセッション cookieと、QueryKeyです。

>>> webenv = search_results["WebEnv"]
>>> query_key = search_results["QueryKey"]

これらの値をsession_cookieとquery_keyという変数に保持するすることで、Bio.Entrez.efetch()のパラメータとしてGIナンバーの代わりにこれらの変数をIDとして使うことができるようになります。

小規模な検索ならば、全てを一度にダウンロードすることに問題はないかもしれませんが、バッチダウンロードしたほうがベターです。restartとretmaxパラメータを検索結果のどの範囲からダウンロードするべきか指定するのに使用することができます。(0からカウントした開始位置と返り値の最大数です。)時にはEntrezからのデータに内部サーバエラー(intermittent error、HTTPError 5XX)が入る場合もあるかもしれません。try: ... except:節を用いてこれに対処します。例：

from Bio import Entrez
import time
try:
    from urllib.error import HTTPError  # for Python 3
except ImportError:
    from urllib2 import HTTPError  # for Python 2
batch_size = 3
out_handle = open("orchid_rpl16.fasta", "w")
for start in range(0, count, batch_size):
    end = min(count, start+batch_size)
    print("Going to download record %i to %i" % (start+1, end))
    attempt = 1
    while attempt <= 3:
        try:
            fetch_handle = Entrez.efetch(db="nucleotide", rettype="fasta", retmode="text",
                                         retstart=start, retmax=batch_size,
                                         webenv=webenv, query_key=query_key)
        except HTTPError as err:
            if 500 <= err.code <= 599:
                print("Received error from server %s" % err)
                print("Attempt %i of 3" % attempt)
                attempt += 1
                time.sleep(15)
            else:
                raise
    data = fetch_handle.read()
    fetch_handle.close()
    out_handle.write(data)
out_handle.close()

目的を説明するため、ここではFASTAのレコードを３つバッチダウンロードしていますが、ゲノム全体や染色体をダウンロードしたいのではない限り、通常はより大きなバッチサイズになるでしょう。

9.15.2:

もう一つ履歴機能の使用例を見せましょう。昨年パブリッシュされたOpuntiaに関する論文を検索し、MedLineフォーマットでダウンロードします。

from Bio import Entrez
import time
try:
    from urllib.error import HTTPError  # for Python 3
except ImportError:
    from urllib2 import HTTPError  # for Python 2
Entrez.email = "history.user@example.com"
search_results = Entrez.read(Entrez.esearch(db="pubmed",
                                            term="Opuntia[ORGN]",
                                            reldate=365, datetype="pdat",
                                            usehistory="y"))
count = int(search_results["Count"])
print("Found %i results" % count)

batch_size = 10
out_handle = open("recent_orchid_papers.txt", "w")
for start in range(0,count,batch_size):
    end = min(count, start+batch_size)
    print("Going to download record %i to %i" % (start+1, end))
    attempt = 1
    while attempt <= 3:
        try:
            fetch_handle = Entrez.efetch(db="pubmed",rettype="medline",
                                         retmode="text",retstart=start,
                                         retmax=batch_size,
                                         webenv=search_results["WebEnv"],
                                         query_key=search_results["QueryKey"])
        except HTTPError as err:
            if 500 <= err.code <= 599:
                print("Received error from server %s" % err)
                print("Attempt %i of 3" % attempt)
                attempt += 1
                time.sleep(15)
            else:
                raise
    data = fetch_handle.read()
    fetch_handle.close()
    out_handle.write(data)
out_handle.close()

執筆時点で、これは２８のマッチを返します。しかしこれは日付に依存する検索なので、当然ながら変化する可能性があります。9.12.1節で詳しく見ています。Bio.Medlineを用いて保存したレコードをパースすることができます。

9.15.3 サイテーションの検索

9.7節でELinkを使用して特定の論文のサイテーションを検索できることを見ました。残念ながら、これはPubMed Centralのためのためにインデックス化されたジャーナルのみをカバーしています。(PubMedの全ジャーナルをカバーするにはNIHの負担がかなり増えます。)ではBiopythonのPDBパーサでPubMed ID14630660をパースしてみましょう。

>>> from Bio import Entrez
>>> Entrez.email = "A.N.Other@example.com"
>>> pmid = "14630660"
>>> results = Entrez.read(Entrez.elink(dbfrom="pubmed", db="pmc",
...                                    LinkName="pubmed_pmc_refs", id=pmid))
>>> pmc_ids = [link["Id"] for link in results[0]["LinkSetDb"][0]["Link"]]
>>> pmc_ids
['2744707', '2705363', '2682512', ..., '1190160']

よし、１１の論文が見つかりました。しかしなぜBiopythonのアプリケーションノート(PubMed ID 19304878)が見つからないのでしょう？変数名から気づいたかもしれませんが、これらは実際のPubMed IDではなくPubMed CentralのIDなのです。我々のアプリケーションノートはリストの３つ目 2682512です。

ではもし、(私のように)PubMed IDのリストが欲しい場合はどうすればよいでしょう？再びELinkを呼び出して、変換することができます。これは２段階のプロセスで行われます。ですので先ほど(9.15節で)説明した履歴機能を使用できます。

しかしまずは、２つ目の呼び出しを(１つ目とは別に)行う、より率直なアプローチを取ってみましょう。

>>> results2 = Entrez.read(Entrez.elink(dbfrom="pmc", db="pubmed", LinkName="pmc_pubmed",
...                                     id=",".join(pmc_ids)))
>>> pubmed_ids = [link["Id"] for link in results2[0]["LinkSetDb"][0]["Link"]]
>>> pubmed_ids
['19698094', '19450287', '19304878', ..., '15985178']

今回は、すぐにBiopythonのアプリケーションノート(PubMed ID 19304878)が３つ目にすぐ見つかりました。

では、ぜひ同じことを履歴機能を用いて行ってみてください。

最後にもう一度、Entrezの呼び出しを行う際はemailアドレスの設定を忘れないこと。

翻訳ここまで

思ったより長くて大変でした…