localBLASTで遺伝子対応表作成
本日(2013年7月25日)、PLOS ONEから論文を発表しました。お世辞にもいい機能アノテーションがなされているとは言えない「非モデル」生物種カイコにおいてマイクロアレイデータ解析を可能とし、ヒトで使われているパスウェイ解析ソフトが利用可能となったことが解析の「突破口」となりました。それを可能としたのはカイコとヒトの「遺伝子対応表」です。
もちろん、あらかじめ個別のモデル生物データベースセンター的なサイトで用意されてあればそれを使うのがいいのですが、必ずしもそれがあれば利用者側で何もやらなくていいというわけではありません。今回の場合もそうだったのですが、そういった提供されているデータセットと、マイクロアレイ解析の場合など調べたいデータセットとのIDの対応関係が簡単にとれない状況であることもあります。そういった場合、やりたい解析(今回の場合、カイコをヒトと比べてこの二種類の生物種間で対応表が作成したかった)に合わせて「対応表」を個別に作成する必要があります。それをどうやったかを手順を追って解説したいと思います。
まず、どこのデータセットを比較対象とするかの取捨選択から解析は始まります。世の中にはすでに利用可能な、たくさんのデータセットがあります。そのうち、どれを使うべきか?これには実は答えがなくて、それぞれの状況や好みに応じて選ぶしかないと思います。UCSCのサイトのものを使うやり方もありますが、Ensemblの方が整然としていて使いやすいので、今回もEnsemblの方を使いました。ここにあるデータダウンロードindexページから必要なデータをダウンロードしてきます。今回の場合、カイコアレイの元のEST(cDNA)配列に対して、ヒトのcDNA配列の対応を付けたかったので、それに相当するHuman cDNA(FASTA)をクリックします。
ftp://ftp.ensembl.org/pub/current_fasta/homo_sapiens/cdna/
クリックすると上記のサイトが開きますが、そこにあるうちallの方、現時点(2013/7/25現在)ではHomo_sapiens.GRCh37.72.cdna.all.fa.gzをダウンロードします。 もちろんFTPサイトのディレクトリ構成を知っているのなら、FTPサイトからftpで直接ダウンロードしてくるのでも構いません。
そして次は、メインの対応表作成の(localな)BLAST実行ですが、localBLASTのインストールは前にも話題にしたのでそちらを参考にしてください。今回の場合の実行コマンドは以下のようになります。
gunzip Homo_sapiens.GRCh37.72.cdna.all.fa.gz
makeblastdb -in Homo_sapiens.GRCh37.72.cdna.all.fa -dbtype nucl -hash_index
tblastx -query kaiko.fa -db Homo_sapiens.GRCh37.72.cdna.all.fa -outfmt 6 -max_target_seqs 1 -evalue 1e-10 -num_threads 4 -out kaiko_human.txt
1行目で、ダウンロードしてきたgzip圧縮されたファイルを解凍します。 2行目で、localBLASTを実行するためのデータベース作成で、解凍したファイルを-inで指定し、-dbtypeには核酸配列なのでnuclを指定します。 3行目がメインのtblastxの実行で、queryにカイコ配列の(FASTA形式の)ファイル(-query kaiko.fa)、dbにダウンロードしてきたヒト配列のファイル(-db Homo_sapiens.GRCh37.72.cdna.all.fa)を指定します。それ以外のオプションも重要で、-out kaiko_human.txtは出力結果のファイル名、-outfmt 6は出力をタブ区切りテキストファイル、-max_target_seqs 1はヒットの最大値が1(つまりトップヒットのみを抽出)、-evalue 1e-10はE-valueのカットオフが1e-10、-num_threads 4はこのプログラムにおいて4CPUまで占有してよい、という指示をコンピュータに与えています。
ここのコマンドを実行する際によく訊かれるのが「E-valueのカットオフ値をいくつにすればいいか?」ということです。これには実は答えがなくて、それぞれの状況に応じて選ぶしかありません。E-valueは「期待値」ですから、1以下であれば統計的には有意なのですが、それが生物学的にも有意かどうかはまた別の問題です。今回、論文のcoauthorにもなっている学生さんといくつかの値をふって条件検討して、出てきた結果が一番生物学的にもっともなものを選んで、1e-10(0.0000000001のことです)としてあります。
すべての計算が終わるまでに時間はかかりますが、この結果カイコのそれぞれの配列に対してヒトcDNA配列の中で(アミノ酸配列に翻訳しながら)配列類似性検索を行なって一番近かったもののリストがkaiko_human.txtというファイル名で出来上がります。
最後に、この作成した対応表を元に、パスウェイ解析ソフト(我々の場合、KeyMolnetを利用しました)でカイコでのマイクロアレイの値をあたかもヒト遺伝子の発現値のように扱うことでパスウェイを描きました。その結果、論文の図にあるようなパスウェイが予測され、その情報を元に定量的PCRで確認の実験をしたわけです。このようなウェットとドライの「キャッチボール」からホームランとはいかないまでもクリーンヒットな成果をこの度出すことが出来て大変うれしいです。今後も頑張っていきたいと思います。
オマケ: FASTA形式ファイルからIDだけを抜き出す簡単なワンライナーを紹介しておきます。解凍したFASTA形式のファイルを使って、
perl -ne 'print "$1n" if(/^>(w+)/)' Homo_sapiens.GRCh37.72.cdna.all.fa > ENST_all.txt
のようなperlスクリプトでFASTA形式ファイルのタイトル行にあるIDを一覧表(ファイル名:ENST_all.txt)として生成することができます。このようなリストはどのIDのものがヒットがあったか(なかったか)を調べるときに必要となるものです。そのうちそのやり方もブログに書きたいと思います。