GGRNAを試す
ある短いオリゴの配列が1000本くらいあって、ソレがどの遺伝子に由来する配列なのかを調べる必要ができた。なんで由来が分からないんだよ！ってのは、まあ事情があるので、しょうがない。
普通に考えるとblast+を使えば良いかと思うんだけれど、フルマッチの配列のみを取り出したいし、GeneSymbolとEntrezのGeneIDを取り出したいという要求もあるので、普通にNCBIからrefseq_rna等をダウンロードして使うのも問題がある*1。
せっかくなので、GGRNAを利用してみることにした。
RESTとJSONを使う
GGRNAは、APIが公開されているので、RESTでrequestして、結果をJSONで得ることができる。
例えば、「GGCGCTAAAAGTTTTGAGCTTCTCAAAAGTCTAGAGCCACCGTCCAGGGAGCAGGTAGCTGCTGGGCTCC」という配列をGGRNAで検索してJSONを得るためには、次のようなrequestを投げれば良い。
perlだとこんな感じ。
my $seq = "GGCGCTAAAAGTTTTGAGCTTCTCAAAAGTCTAGAGCCACCGTCCAGGGAGCAGGTAGCTGCTGGGCTCC"; my $url="http://ggrna.dbcls.jp/api/hs/seq:" . $seq . ".json"; my $browser = WWW::Mechanize->new; $browser->get($url); my $content = $browser->content();
すると、だいたい、以下のようなJSONのテキストが返ってくる。
{"summary":[{"keyword":"seq:GGCGCTAAAAGTTTTGAGCTTCTCAAAAGTCTAGAGCCACCGTCCAGGGAGCAGGTAGCTGCTGGGCTCC","count":"5"},{"keyword":"[AND]","count":"5"}],"version":"GGRNA : RefSeq release 54 (20120713)","time":"2012-09-03 15:59:43","results":[{"snip":"..... cctcccatgtgctcaagactggcgctaaaagttttgagcttctcaaaagtctagagccaccgtccagggagcaggtagctgctgggctccggggacactttgcgttcggg .....","geneid":"7157","symbol":"TP53","version":"NM_000546.5","sequence_position":"position 40 (CDS: 203 - 1384) ","definition":"Homo sapiens tumor protein p53 (TP53), transcript variant 1, mRNA.","refid":"NM_000546"},{"snip":"..... cctcccatgtgctcaagactggcgctaaaagttttgagcttctcaaaagtctagagccaccgtccagggagcaggtagctgctgggctccggggacactttgcgttcggg .....","geneid":"7157","symbol":"TP53","version":"NM_001126112.2","sequence_position":"position 40 (CDS: 200 - 1381) ","definition":"Homo sapiens tumor protein p53 (TP53), transcript variant 2, mRNA.","refid":"NM_001126112"},{"snip":"..... cctcccatgtgctcaagactggcgctaaaagttttgagcttctcaaaagtctagagccaccgtccagggagcaggtagctgctgggctccggggacactttgcgttcggg .....","geneid":"7157","symbol":"TP53","version":"NM_001126113.2","sequence_position":"position 40 (CDS: 203 - 1243) ","definition":"Homo sapiens tumor protein p53 (TP53), transcript variant 4, mRNA.","refid":"NM_001126113"},{"snip":"..... cctcccatgtgctcaagactggcgctaaaagttttgagcttctcaaaagtctagagccaccgtccagggagcaggtagctgctgggctccggggacactttgcgttcggg .....","geneid":"7157","symbol":"TP53","version":"NM_001126114.2","sequence_position":"position 40 (CDS: 203 - 1228) ","definition":"Homo sapiens tumor protein p53 (TP53), transcript variant 3, mRNA.","refid":"NM_001126114"},{"snip":"..... cctcccatgtgctcaagactggcgctaaaagttttgagcttctcaaaagtctagagccaccgtccagggagcaggtagctgctgggctccggggacactttgcgttcggg .....","geneid":"7157","symbol":"TP53","version":"NM_001126118.1","sequence_position":"position 40 (CDS: 437 - 1501) ","definition":"Homo sapiens tumor protein p53 (TP53), transcript variant 8, mRNA.","refid":"NM_001126118"}]}
このJSONを正規表現を使いこなしてparseする必要は...ない。CPANにJSON関連のパッケージが大量にあるから、それを使えばいいはずだ。だけれど、何を使ったら良いのかよく分からない。雑に「use JSON;」だけしてみたり、こっちのほうが早いという評判を聞いて「use JSON::XS;」で実行してみると、何故かmalformed JSONと呼ばれる。よく分からないので、検索した結果からコピペしてみた*2ら、ちゃんとparseできた。
use JSON -support_by_pp; my $json = JSON->new; my $json_text = $json->allow_nonref->utf8->relaxed->escape_slash->loose->allow_singlequote->allow_barekey->decode($content); my %genes=(); foreach my $c(@{$json_text->{results}}){ $genes{$c->{geneid}}=$c->{symbol}; } while(my ($key, $value) = each (%genes)){ print "$key = $value\n"; }
コレを実行すると、上記の配列のGeneSymbolがTP53でそのGeneIDは7157だということが分かる。
実行してみた結果
レスポンスも一瞬ですばらしい。これはコピペスクリプトが良いのもあるんだろうけれど、何よりもサーバーのレスポンスがイイんだろうと思う。自前のダサいPCでstandalone blastを動かせるようにしてああだこうだするよりも、ずっとイイ。もっと前から知っておくべきだった。
んで
このGeneIDを使って、あらかじめダウンロードしてsqlite3のDBにつっこんでおいたgene_infoやらgene2accessionやらgene2goやらgene2ensemblで芋づる式にあれこれデータを引っ張れば、あっという間に「ソレっぽいアノテーションテーブル」の出来上がり。
