2011年6月2日木曜日

発現解析パイプラインを作るぞ！　その1：　TopHat　の使い方

TopHatといえば、NGS発現解析で良く使われるツールです。
論文やポスターでも頻繁に見かけるので、一応、スタンダードなツールと見て良いでしょう。
TopHatと一緒にCufflinksという名前も、聞いたことがあるのではないでしょうか？
こちら、セットで使うと、RNA-Seqの解析にそこそこ威力を発揮します。

ちなみにこの写真は、エイブラハム・リンカーンが、実際に使っていたTopHatです。
スミソニアンのアメリカ史博物館での一枚です。
残念ながらCufflinks（カフスボタン）は、飾っていませんでした。
冗談はさておき、TopHatは、Exon-Intronのジャンクションを挟んでうまくアライメントするマッピングツールです。　Cufflinksは、そのマップされたフラグメント（リード）をカウントするツールです。
Cufflinksはこちらから（http://cufflinks.cbcb.umd.edu/）
ダウンロードできます。　マニュアルもあります。　Cufflinksは最近頻繁にバージョンアップされています。
TopHatを使うには、Bowtie（高速マッピングツール、http://bowtie-bio.sourceforge.net/index.shtml）
と、Samtools （SAM<->BAMファイルの変換やSNP検出などをするツール、http://samtools.sourceforge.net/）も一緒にインストールしておくと便利です。

Bowtie、TopHat、Samtoolsは、インストールも比較的難なくできますが、Cufflinksはちょっと準備が必要（bjamとかいうエンジンを入れる必要）ですので、行き詰ったらLinuxに詳しい方に助けを乞いましょう。　

4つともうまく入ったら、どんな風に動くのか、早速試してみましょう。
デモデータは、SRAから落としてきたSRP003186のデータを使います。　乳がんサンプル、イルミナGAIIのペアエンドデータです。
そうそう、SRAのファイルは、そのままでは解凍できないので、SRAのサイトからSRA-TOOLSというものを落としてきて、Linuxマシンにインストールします。　詳しくは昨年12月の書き込み（サンプルデータの取得法　２　NCBI 、http://shortreadbrothers.blogspot.com/2010/12/ncbi.html）にあるが、ここでのコマンドはfastq-dump2 から現在 fastq-dump に変わっているので注意！


全体の流れを示します。既知トランスクリプトーム情報を使った、発現変動解析パイプラインの大まかな例です。　必要以外のジョブもありますが、そこはつっこまないで。
  1. リファレンス配列Hg19を用意する
  2. Hg19に対してBowtieのインデックスを貼る
  3. TopHatを実行、リファレンスにマッピングする
  4. できたBAMファイルをソートしておく
  5. 後でViewerで見やすいようにソート後のBAMファイルにインデックスを貼る
  6. ついでにソートしたBAMファイルをSAMファイルにもしておく
  7. そのまたついでにソートしたBAMファイルをもとに、リファレンス配列と比較して異なる塩基（変異）のポジションを見つける
  8. さっき4でソートしたBAMファイルから、既知トランスクリプトのGTFファイルをもとに発現比較する ― Cuffdiff 実行
さて、リファレンス配列は、UCSCのFTPから適当なfa.gzをダウンロードしてきて、展開。
24本の染色体だったら、24個のファイルができるので、それを　cat *.fa > hg19.fa としてひとつのファイルにまとめる。
ついでに、UCSCから、好きなアノテーションファイル（GTFフォーマット）を落としてきます。
これについては、この画面を参照　

  
インデックスは、Bowtieのbuildコマンドで
bowtie-build /Path_to_reference/UCSC_hg19/hg19.fa ./bowtie_ref/hg19

いよいよTopHatの実行です。
先ほどのgtfファイルを用意して、以下のようなコマンドを流します。　3つのサンプルがあります。
コマンドは3行です。ペアエンドなので、それぞれ***_1.fastq ***_2.fastq　というファイル名です。

tophat -r 100 -p 4 -G /Path_to_gtf/refGene.gtf -o ./tophat_SRP003186_Control /Path_to_bowtie_ref/hg19 /Path_to_read_file/Control/SRR064437_1.fastq /Path_to_read_file/Control/SRR064437_2.fastq
tophat -r 0 -p 4 -G /Path_to_gtf/refGene.gtf -o ./tophat_SRP003186_MCF7/Path_to_bowtie_ref/hg19 /Path_to_read_file/MCF_7/SRR064286_1.fastq /Path_to_read_file/MCF_7/SRR064286_2.fastq
tophat -r 100 -p 4 -G /Path_to_gtf/refGene.gtf -o ./tophat_SRP003186_SKBR3 /Path_to_bowtie_ref/hg19 /Path_to_read_file/SK_BR_3/SRR064441_1.fastq /Path_to_read_file/SK_BR_3/SRR064441_2.fastq

パラメータについて
-r: 　inner mate distance 平均フラグメントの長さ－両側のリードの長さ。MCFの例の場合、フラグメント長は100塩基でリード長は50塩基なので、100-50x2=0
-p：　使用するCUP数（デフォルト：1）
-o：　アウトプットファイルを書き出すディレクトリ
-G：　gtfファイルの場所
このほかのパラメータについてはおいおい。　もちろんTopHatのウェブサイトにも情報はたくさんあります。

アウトプットファイルは、たくさんできるのですが、以下の3ファイルがとくに重要！というか、後はあまり見ません。
  1. accepted_hits.bam：BAM フォーマットのアライメントファイル
  2. junctions.bed：Exon Junction情報をまとめたファイルで、UCSC BED track フォーマット。ジャンクションは2つのBEDブロックから成り、"maximal overhang" で指定した数以上の塩基によって構成しているはず。スコアはアライメントの数を示す。UCSC ゲノムブラウザに表示できる
  3. insertions.bed / deletions.bed：insertions/deletions情報をまとめたUCSC BED tracksファイル。
    Insertions.bed の chromLeft は、Insertion直前の塩基を表し、
    Deletions.bed の chromLeft は、Deletionする最初の塩基を表している、らしい
 さて、ここまでいかがでしょうか？

つづく
投稿者 時刻:
ラベル: ,

0 件のコメント:

コメントを投稿

登録: コメントの投稿 (Atom)