配列セットを取得する(TogoWS編)
cDNA配列のセット(多くの場合数万のオーダーの配列群)が論文発表されたのに、それがウェブサイト上にないとかいう話をよく聞きます。その場合、大変申し訳無いのですが論文の読み込みが足りないと言わざるを得ません。論文発表されたからには必ずデータは公共データベースに公開されているはずです。していないのは「公開した」とは言えないからです。 ＃まれにNGS黎明期のゴタゴタで配列を公開していないのに論文が通ってしまっている例がありますが、それは例外として。
論文にはデータベースのアクセッション番号がどこかに書いてあるはずです。まずそれを探しましょう。例えば、
All full-length cDNA sequences obtained in this work appear in public databases under accession nos. [DDBJ:AK377185-AK388575].
とあると、アクセッション番号AK377185-AK388575で公開しているわけですから、以下の様なPerlのコード(get.prl)で配列が取得できます。 [perl] my $togows = "http://togows.dbcls.jp/entry/nucest"; #baseURL my $start = 377185; my $end = 388575; foreach my $i ($start..$end) { my $id = "AK$i"; #accession No. print STDERR "getting $id...n"; system("curl $togows/AK$i.fasta 2> /dev/null"); sleep 5; } [/perl] 2,3行目で取得開始と取得終了の番号を指定しているほか、4-9行目でその間配列取得を繰り返す、そんなものすごくシンプルなコードです。 7行目が実際に配列を取得する部分で、システムにcurlコマンドがインストールされていることが前提となっています(MacOSXには最初から入っています)。DBCLS謹製のTogoWSを使うと、 [html] http://togows.dbcls.jp/entry/nucest/AK377185.fasta [/html] の様なURLを叩くとAK377185の核酸配列がfasta形式で取得できます。それを利用して、各データベースエントリの取得を繰り返す、ただそれだけです。
TogoWSでどのようなデータベースや形式をサポートしているかはREST サービスのサンプルを御覧ください。
連続してい一度に取得するとサーバー側に負担をかけますので8行目にあるように1エントリ取得するごとに5秒sleepして繰り返すようになっています。
そして、 [shell] perl get.prl > fasta.txt [/shell] とすればfasta.txtに欲しい配列が約5秒毎に追加書込されていきます。
これの実行はそれなりに時間がかかります。なぜなら1エントリ取ると5秒休みがあるわけですから、配列取得が一瞬で終わったとしても1分で12エントリ、1時間で720エントリほどしか取得できないからです。そこで途中どれぐらい取れたか、が気になると思うのですが、出力形式がFASTA形式ですから、 [shell] grep -c ^> fasta.txt [/shell] とすると先頭が>で始まるヘッダー行の数を数えることができ、どれぐらいの配列が取得できているかがわかります。
繰り返しますが、大事なのは論文などの文書中に書かれたアクセッション番号のリストを読み解いてコンピュータにそれを伝えることです。それが人間のやる仕事です。同じ操作を一万回以上も繰り返してやるのはコンピュータの仕事です。お互いうまく住み分けて効率的に研究を進めていきましょう。