Wolf Ears

←前回　Tophat・Cufflinksを用いた遺伝子発現解析の方法 (3)

前回ではTophatを用いて次世代シーケンサから得られた配列をマッピングしましたが，今回はCufflinksを用いて発現量の推定を行ないます．ひとまず，ここまでやってきたことに加えてこれからやるCufflinksで処理するという手順が，一連の遺伝子発現解析における流れだと思います．もちろん目的に応じて方法は変えていく必要はありますが，人が見ても理解できないような生のデータを生物学的に意味のあり議論の題材となるような標準形式に変換するという点では，ひとつの重要な作業だと思います．このあたりが一番機械的な作業であり，また繊細なチューニングを要するところなのでしょうが，取り敢えずはさらりと表面を流す程度で行きたいと思います．

まずはCufflinksの準備です．CufflinksもTophatと同様にサイトから入手します．
http://cufflinks.cbcb.umd.edu/tutorial.html
通常はReleaseのところにあるpre-compiled binaryを使用してください．もしソースコードからコンパイルする場合は，C++のライブラリであるboostのパスを通す必要があるので少しメンドくさいです．
なお，このCufflinksにはcufflinks,，cuffcompare，cuffdiffといった実行ファイルが入っています．これ以降述べるcufflinksであったりcuffcompareはCufflinks内の個別の実行ファイルのことを指すものとします．

Cufflinks
それでは，Cufflinksを使って発現量の推定を行ないます．

cufflinks -p 2 accepted_hits.sam

発現量推定に必要なマッピングの情報はすべてaccepted_hits.samファイルに入っているので，何かファイルを追加したりして処理する必要はありません．ただ単純にTophatから出てきたものをcufflinksに流すだけです．ただ，細かな設定として実験環境に合わせてオプションを設定する必要があります．まず，上のコマンドに書かれている-pはTophatと同様に使用するスレッドを指定します．他にも代表的なオプションとして，-mはペアエンドで実験しているときのexpected (mean) inner distanceを指定するといったものがあります．
さて，cufflinksをかけると以下のファイルが出力されます．

• transcripts.gtf
• transcripts.expr
• genes.expr

ここで重要なファイルはtranscripts.gtfです．このファイルに，実際に染色体のどの位置でどれくらい発現しているのかといった情報が含まれています．このファイルは拡張子の通りGTFと呼ばれるフォーマットとなっており，GFFと同様に遺伝子のアノテーション情報を記載するときの標準的なフォーマットとなっているようです．
他にもtranscripts.exprなどのファイルも出力されますが，これらのファイルはCufflinksのmanualによると

these are reserved for debugging, and may change or disappear in the future

とあるので，あまり気にする必要はないかと思います．

Cuffcompare
それでは次にCuffcompareを使って，今回得られた発現領域や発現量と既存のReference Genome sequence annotation gtfとの比較を行ないます．

cuffcompare -r reference.gtf transcripts.gtf

ここで，必要になるのが比較する元となるアノテーションファイルです．これはマイクロアレイなどの実験などによって求められた発現位置が遺伝子名などとともにタグ付けされたものです．このアノテーションファイルを今回得られた領域と比較して対応付けることで，どのようなトランスクリプトームが発現しているかを知ることができます．
まずはこのアノテーションファイルを入手します．ファイルはNCBIのものとEnsemblのものの2種類がありますが，大きな違いとしては遺伝子IDなどの表記の差なので，基本的にはどちらでも構わないと思います．もし，どちらかのアノテーションで統一する必要があったり，個人的に使いやすい方があればそちらを選べばいいと思います．ただ，どれくらいキチンとアノテーションされているかが少し違うようです（またいずれ記事を書きたいです…）．
Ensemblのアノテーションファイルについては，大抵はftpサイトなどに置かれているようです．また，NCBIのアノテーションファイルについては，UCSCのTable BrowserからRefseq Genesというtrackを指定することでダウンロードすることが出来ます．
ftp://ftp.ensembl.org/pub/current_gtf/
http://genome.ucsc.edu/cgi-bin/hgTables
一応まとめっぽいサイトも載せておきます．
http://wiki.geneontology.org/index.php/Reference_Genome_sequence_annotation

これでアノテーションファイルの準備は完了です．あとはCuffcompareをかけるだけです．
Cuffcompareをかけると以下のファイルが出力されます．

• transcripts.refmap
• transcripts.tmap
• その他諸々

transcripts.refmapは，リファレンスのアノテーションに対応するtranscriptsを表します．言い換えると，リファレンスのtranscriptsがどのCufflinksで推定されたtranscripts対応しているかということを表しています．一方，transcripts.tmapは，Cufflinksから推定されたtranscriptsに対応するリファレンスのtranscriptsを表します．つまり，これが今回推定したトランスクリプトームとリファレンスのトランスクリプトームの一対一に対応付けられたファイルとなります．
少々ややこしく分かりにくい説明となってしまいましたが，要するに

reference transcripts <-> cufflinks transcripts

のどちらを基準として対応付けしているか，ということです．

さて，ここで関心があるのは，自分が出した結果がリファレンスのアノテーションとどのように紐付けされているか，つまりtranscripts.tmapの方のファイルです．transcripts.tmapの中身を覗いてみると，cuff_idの順番で一行ごとにリファレンスのアノテーションと対応付けされたトランスクリプトームがタブ区切りのリストとなっています．ここで，ref_gene_idやref_idがリファレンスのアノテーション由来のものとなっており，”-”となっている箇所は該当するものがなかったということを表しています．次に，ここで一番重要となるのがFPKMという値です．これは，Fragments Per Kilobase of exon per Million fragmentsの略で，RPKMなどと同様に発現量の一種の指標となります．
では，試しにこのFPKMを基準ににソートして，発現量が多い順に並べ替えてみます．

sort -k 7 -nr transcripts.tmap > output.tmap

これでFPKMの値が大きい順にソート出来ます．あとは煮るなり焼くなり…


といった感じで，取り敢えず「Tophat・Cufflinksを用いた遺伝子発現解析の方法」は終了です．これまで駆け足で各ステップの流れを大まかに説明してきましたが，かなりまとまりのない雑多な内容になってしまった感じはあります．それも，このシリーズの元となる部分はここ2，3ヶ月の間やってきたことが元となっているのですが，あまりにも事前知識なく手探りな状態でここまで進んできたことが主な原因ですね．かなりアバウトさが滲み出る結果となりました．ただ，これを書くにあたって，これまで曖昧にしてきた部分などが浮き彫りになり自分の勉強にもなったので，その点は良かったと思います．これからも，細かなことや間違っている部分は訂正し，足りないと感じることがあればこれらの記事にも直接書き足していく予定ですので，よろしくお願いします．