本業の方で色々忙しくなっておりまして、更新の間が開いてしまいました。
今回はRNA-seqについて語りたいと思います。 RNA-seqはChIP-seqよりもメジャーなので、日本語での解説ブログも充実していますが、情報が古いものだと今だにtophat-cufflinksを使っていたりします。それじゃだめだよ！というお話。
ここではツールのレビューがメインです。RNA-seq解析のHowtoそのものについては下記サイトが大変わかりやすいですので、参照されてください。RベースでのRNA-seq解析法です。
RNA-seq
RNA-seqの目的は大きく分けて、
- 既知遺伝子の発現量を網羅的に計測し、複数サンプルで発現変動比較する
- 未知の転写物（isoform含む）を同定する
の２つがあり、前者のケースが大半かと思います。
前者の場合、既知RNAのデータベースに対してリードをマッピングします。後者の場合はゲノムに直接マッピングした後アセンブリングするか、まずリードをアセンブリングしてからゲノムにマッピングすることになります。
RNA-seqでの発現量解析
の３ステップに分けられ、それぞれ異なるプログラムが必要です。 ステップの詳細が知りたい方はConesaらの論文*1を参照されると良いでしょう。
Tophat-cufflinks
Tophatはリードをマッピングするツール、cufflinksは発現量計測、発現変動解析、新規転写物の同定などのコマンドがセットになったプログラム集です。
tophat-cufflinksは様々な目的に利用可能であること、解説ページが豊富に公開されていたこと、1 replicateでも発現変動比較ができるといったメリットがあり、一時期はかなり普及していました。
その後tophat-cufflinksは開発・サポートが終了し、後継のhisat2またはkallistoに切り替えることが開発者より推奨されています。精度面の問題も明らかになってきており、今ではほとんど使われることもなくなったのですが、論文の査読をしているとtophat-cufflinksを使っているものに未だに出くわします。初心者ならともかく、ある程度サーベイして知識を持った人が未だに使っているのは正直勉強不足だろうと思いますが。。。
最近のツール事情
RNA-seqのツール開発は今でも群雄割拠の様相を呈しており、ツールごとの精度比較を行ったサーベイ論文がいくつか発表されています*2*3。この論文によれば、cufflinksはsingle exonの遺伝子に対して非常に精度が低くなることが報告されています。また、single-end RNA-seqでisoformレベルの発現量を推定した場合にも精度が悪いようです。また、１サンプルからの発現変動比較が可能ですが、１サンプルだと当然精度が悪いです。 従って、paired-endを使って遺伝子レベルの発現量比較をしている場合のみ、tophat-cufflinksはいちおう今でも利用可能ですが、特に使い続けるメリットはありません…査読者にも必ず突っ込まれるでしょう。
ツールの推奨組み合わせ
基本的にはマッピングー発現量計測ー発現変動解析の流れで用いるツール群はセットになっています。 検証されていない組み合わせを用いると予期せぬエラーが起きる可能性がありますので、素直に推奨組み合わせを使いましょう。 具体的には以下のような組み合わせがあります（/ で挟んでいるものはどちらも利用可という意味です）。
- (bowtie2/STAR)ーRSEMー(edgeR/DESeq2)
- Hisat2ーStringTieーBallgown
- kallistoーsleuth
- Salmonー(edgeR/DESeq2)
いわゆる正統派の「マッピングー発現量計測ー発現変動解析」という意味では、 (bowtie2/STAR)-RSEM-(edgeR/DESeq2)を使うのが良いと個人的には思います。精度は最高レベルで、single-endでも利用可、遺伝子単位にもtranscript単位にも対応しており、利用実績も多いです。bowtie2とSTARは精度はそれほど変わりませんが、STARの方が高速です（そのかわりメモリ消費がとても大きい）。edgeRとDESeq2はどちらも優れたツールなので、お好みでよいと思います。他にもいくつかツールはありますが、精度面でこれらを大きく超えることはありません。TIGAR2という国産の発現量計測ツールもありますが、ものすごく遅い割に精度はそんなに変わりません。。
Hisat2-StringTie-Ballgownは新規転写物の同定が可能で、プロトロルがNature protocolsに公開されています*4。一方、既知の遺伝子発現計測では精度はそれほどでもないです。 グラフゲノムアラインメントというコンセプトが取られており、最近はメタゲノム（細菌叢）の解析への応用が検討されています。 なお、eXpressという類似のツールがありますが、こちらはサポート終了していますので使わないようにしましょう。
kallistoとSalmonは"pseudo-alignment (偽アラインメント)"と言って、リードを遺伝子にひとつひとつマップすることなく遺伝子発現を計測可能となったツールです。そのため大変高速です（試してみるとびっくりします）。これはひとつのエポックメイキングであり、精度面でも最高レベルに近いため、「これからはpseudo-alignment だろう！」という機運はあるのですが、デフォルトではBAMファイルを生成してくれない（マッピングしないので）ことや、利用実績が少なくちょっと怖いこともあり、今のところ様子見している人が多いように感じます。なおSailfishというツールもありますが、これはSalmonの先行ツールです。
kallistoは少し癖があります。基本的にはpaired-endのみに対応で、sleuthという独自の発現変動解析ツールがペアになっており、これ以外は使わないことが推奨されています（組み合わせた場合にどうなるのかは未確認です）。そのためedgeRなどを使いたい人は使えないということになってしまいますね。
また、kallisto・Salmonともtranscript単位でしか発現量を出力してくれません。遺伝子単位での発現量マトリクスを生成したい場合はtximport*5というツールを併せて使いましょう。
まとめ
ここではマッピングー発現量解析ー発現変動解析のためのツールを簡単に紹介しました。 他にも、正規化手法やread assembling, alternative splicingなどさまざまな目的の手法が本当にたくさんあり、全てサーベイするのは難しいほどです。
ちなみに私はSTAR-RSEM-edgeRを使っています。pseudo-alignmentの高速性は大変魅力的なのですが、万一後から大問題が発覚すると怖いので、今のところ保守的に行っとこうという感じですね。。。お試しレベルではkallistoもSalmonも使ってはいます。おそらく二年後くらいにはまた様子がすっかり変わっているのでしょうね。
*1:A. Conesa et al., A survey of best practices for RNA-seq data analysis, Genome biology 2016, DOI 10.1186/s13059-016-0881-8
*2:A. Kanitz et al., Comparative assessment of methods for the computational inference of transcript isoform abundance from RNA-seq data, Genome Biology, 2015, DOI 10.1186/s13059-015-0702-5
*3:M. Teng et al., A benchmark for RNA-seq quantification pipelines, Genome Biology, 2016, DOI 10.1186/s13059-016-0940-1
*4:M. Pertea et al., Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown, Nature Protocols, 2016, doi:10.1038/nprot.2016.095
*5:http://bioconductor.org/packages/release/bioc/html/tximport.html
