このチュートリアルは、生物学・医学・農学系研究者がRNA-Seqのデータ解析を学習することを想定し、発現量データ（数値のテキストファイル）のインポートから生物学的解釈まで一通りの流れを解説したものです。FASTQファイルから始める場合、FASTQから発現量を算出までのパイプラインを別のチュートリアルにまとめてありますのでそちらをご覧ください。
ここではSubio Platform という解析ソフトを使っています。Rに比べて知名度はありませんが、学会発表や論文投稿にももちろん使えます。初心者がR/Bioconductorで学ぼうとすると、コマンドの習得に集中しすぎて、肝心のデータ解析に対する理解が進まないことがよくあります。Subio Platformでは、その操作がデータにどのような効果を与えるか、どのような特徴を抽出するかを一つ一つ視覚的に確認しながら進めることができるので、データ解析をきちんと学びたい初心者に最適です。たとえRで解析するとしても、一つ一つのコマンドがどのような意味を持つのかを理解していることはとても大事です。一度Subio Platformを使ってこのチュートリアルを実践してみるのは有益でしょう。
１．RNA-Seqの遺伝子レベルの発現データのインポート
GSE49110 は、８サンプルから成るRNA-Seqデータです。 このデータをダウンロード して解析してみましょう。このデータセットのSSAファイルはこちらからダウンロードできます。ここでは、Countのデータをインポートしていますが、FPKMやTPMのテーブルの場合も同様にインポートできます。
まずはダウンロードしたファイルをエクセルで編集しますが、ここにあるとおりちょっとしたワザが必要です。いらない列と行を削除して、IDである遺伝子名と発現量のテーブルに整形します。
この編集したテキストファイルを Subio Platform にインポートします。Import Samples ウィザードを開始したら、最初のページで “Multiple Samples in One File” オプションを使い、次のページで “Create A New Platform” を使って下さい。
インポートが完了したら、SOFT formatted family file からサンプルの属性情報を取り込んでおきます。これでより多くの情報を見ることができ、また、キーワードでフィルターをかけられるようになります。
それでは、この8個のSampleから成るSeriesを作って、データの視覚化および解析をしましょう。
２．Series の作成と設定
Seriesがロードされたら、画面左のSeriesパネルにMeasurement ListやDataSetなどのオブジェクトが表示されます。Analysis Browserの上段ではデフォルトでScatter Plotが描画されます。一方、下段ではSetup Seriesタブが開きます。
Setup Seriesタブでは、基本的に画面の左から順にボタンを押していけばいいですので、最初のボタンEdit Parametersを押します。通常は、Sample Informationから情報をインポートするのが早いです。このデータセットは、４つの状態（コントロールと３種類のsiRNA処理）があり、それぞれ２回の反復があります。
Setup DataSetタブに移動し、DataSetの編集と作成を行います。先ほど設定したParameterの優先順位を設定することでSample Groupを定義・整理するのがこのタブの役割です。このチュートリアルでは２つのDataSetを設定します。さらに、Sample Info.タブに情報を記録し、関連する添付ファイルを保存します。
３．正規化と前処理
正規化および前処理は、これをどうするかによってその後の解析結果が変わってしまう大事なステップです。しかも、実際のデータは教科書が想定しているようなものとは限りません。目の前のデータの特徴を正しく理解し、それにあった処理をしなければ、その後の解析で誤った結論を導き出してしまうので十分に注意して下さい。
このチュートリアルでは、プリセットの「RNA-Seq (Counts)」というシナリオから出発して、データに合わせて調整します。その過程で、何を知るためにどこを見るか（操作するか）、そしてどのように判断するかの一例を学んでください。繰り返しますが、ここで紹介するやり方がどんなデータにも通用するわけではありません。実際のデータ解析では、そのデータに合わせてその都度調整が必要です。自信を持てない場合は、どうぞデータ解析サービスをご利用ください。
こちらを併せてご覧ください。
４．フィルタリング(Quality Control)
この前の正規化の説明で見たとおり、測定値があったとしてもすべての測定値が信頼できるわけではありません。解析する前に、解析に値する測定値を抽出し、解析不能の値は解析の邪魔にならないように操作する必要があります。ここで使うのが、Basic Plug-inに含まれるFilterツールです。プラグインをお持ちでない方は、5日間の無料お試しをご利用ください。
フィルターを使う際の基本的な考え方は、「信頼できる遺伝子を抽出する」ではなく、「解析に値しない遺伝子を抽出して除く」です。この二つの違いは解析経験がないとまず分からないでしょう。前者は、例えばコントロールでは発現していなかったのに、何らかの処理をすると発現してくる遺伝子を逃してしまうことになります。この違いを、フィルターツールを実際に使ってみて確かめてください。
もう一つの基本的な考え方は、２段階で行うということです。まず、値が低すぎる遺伝子を除きます。次に、発現変動しない遺伝子を除きます。
ところで、RのDESeq2を使った解析手順を独学していると、次のようなサンプルコードがあります。
dds <- dds[rowSums(counts(dds)) >= 10,]
「全サンプルで合計10リード未満の遺伝子をフィルタリングする」というのは適切でしょうか？経験かまったくまったく十分ではないと思いますし、また、このやり方はサンプル数にも依存します。６サンプルを解析する場合と、６０サンプルを解析する場合では、単純計算でも閾値は10倍違うことになりますよね。では適切な閾値とはどうやって決めたらいいのでしょうか？その答えは、データを見てみないと分かりません。だからこのチュートリアルのように、データをみて理解しながらどのような手順で解析するかを決めていく必要があるのです。
５．PCAとクラスタリング
前段までで準備が整いましたので、ここからいよいよ、さまざまな解析手法を駆使してデータから何らかの意味を抽出する段階に入ります。すべての遺伝子を対象とするのでなく、フィルターを使ってQuality Controlを通過した遺伝子群だけを対象とすることにご注意ください。
まずは発現プロファイルの全体像を俯瞰することで、データの概観を把握し、解析のポイントを明らかにするのがよいでしょう。この目的で役に立つのが主成分分析（PCA）です。PCAの結果を見るときの要点は３つです。まず、点と点の距離が近ければ発現プロファイルが似ている、そして離れていれば発現プロファイルも大きく異なることを表します。次に、原点からの方向です。同じ方向であればある遺伝子群が似たような動きをしていることを表し、原点から遠いほど変動幅が大きいことを表します。また、原点を挟んで点があるということは、発現変動が逆向きであることを示しています。そして、最後に主成分と寄与率です。縦方向と横方向は、あるいはそれぞれの主成分は異なる遺伝子群の動きを代表していると言えます。寄与率が大きいほど、そのような動きをする遺伝子の数が多く、概観を捉えていると言っていいでしょう。しかし、寄与率が大きいことが必ずしも生物学的に重要というわけではありません。生物学では、ごく少数の遺伝子群の動きが将来の方向性を決定づけることはよくあります。つまり、寄与率の小さな主成分が決定的な動きを捉えているかもしれません。
以上を踏まえて、このデータのPCAの結果を見てみましょう。繰り返しの２点間の距離は、siRNA間の距離に比べてずっと小さく、ばらつきの小さい良質の実験データであることが分かります。原点のコントロールと比べて、各siRNA処理群は同一方向（右下）への動きと、個別の動き（右と下）があるようです。
次に、階層型クラスタリングを適用してみます。ヒートマップを概観すると、すべてのsiRNAで共通で発現上昇（赤）する遺伝子は多くあるのと対照的に、発現が下降（青）する遺伝子はsiRNAごとに異なっていることが分かります。解析者がすべきことは、発現差のある遺伝子リストや図を作ることではなく、このような違いや特徴に気づき、細胞内で何が起きているかを洞察することです。
また、PCAと階層型クラスタリングは、基本的には同様の解析結果を違う見方をしているにすぎません。二つの結果を見比べながら、データの理解を深めましょう。
６．発現差のある遺伝子の抽出
このケースでは、前章のクラスタリングで見た通り、発現が上昇する遺伝子群は3種のsiRNAで共通なのに対し、発現が下降する遺伝子群の共通項はあまりに少ないという偏りがありそうです。Basic Plug-inのツールを使って発現差のある遺伝子を抽出したら、Venn Diagramツールで確認してみましょう。
７．遺伝子アノテーションと、エンリッチメント解析
RNA-Seqのデータの場合、countまたはFPKMのテーブルにGene Symbol（ENSGやENSTのID、Entrez Gene IDの場合もある）しか付いていないことも多いです。そこで、遺伝子アノテーションのテーブルをデータベースサイトから取得して埋める必要があります。ここではNCBI FTP Siteサイトを使いますが、ENSGやENSTのIDがついている場合は Ensembl BioMartから取得してください。それ以外のIDの場合は、そのIDと遺伝子アノテーションを管理しているデータベースサイトを使ってください。
遺伝子アノテーションのインポートが終わると、Subio Platformでの検索機能や、Advanced Plug-inに含まれるEnrichment Analysisツールが活用できるようになります。
用語が紛らわしいのですが、Gene Ontology (GO) 解析、パスウェイ解析、ネットワーク解析などのキーワードで尋ねてこられる場合、その多くがエンリッチメント解析のことをおっしゃっています。Gene Set Enrichment Analysis（GSEA）、Ingenuity Pathway Analysis (IPA)、David Functional Annotation、Meatscapeなどはエンリッチメント解析を行うツールです。
エンリッチメント解析については、別に詳しいチュートリアルがありますので、併せてご覧ください。
８．ゲノム上の位置特異的に発現制御されている遺伝子と、モチーフ配列
遺伝子アノテーションのうち、染色体上の位置情報は特別に用意された専用の列に、決まった形式で入れてください。こうするだけで、Subio PlatformではGenome ViewとChromosomeタブが使えるようになります。
RNA-Seqのデータは、原理的には遺伝子ごと、あるいは転写物ごとの発現量を推定することができます。しかし、転写物単位での発現量を推定するということは、リード数が転写物ごとに割かれることになるので、一般的なリード数ではダイナミックレンジが狭くなり、データの質が低下します。従って、現実的には遺伝子単位の発現量の推定値のみが解析対象として有用です。遺伝子単位で解析するということは、転写開始位置等を正確に反映した解析はできないということを承知しておく必要があります。
発現差のある遺伝子を抽出したり、クラスタリングで特定の発現パターンの遺伝子を抽出したら、その遺伝子群がゲノム上に偏って存在していないかどうかを確認することができます。もし特定の領域の遺伝子群が一斉に発現上昇または下降していたら、その領域のエピジェネティックな状態変化や、染色体の構造変化がその発現変動の原因になっている可能性が出てきます。特に偏りが無ければ、転写因子によって発現制御されている可能性が高いと考えられます。
一方、ゲノム配列のデータをダウンロードして、gzで圧縮されたままで構わないのでローカルに保存していれば、モチーフ配列の位置を検索できます。モチーフ配列の検索は、IUPAC nucleotide codeに対応しています。モチーフ配列の位置が取得できたら、それを転写開始位置近傍に持つ遺伝子を検索して、それら遺伝子の発現パターンを確認することができます。
このケースでは、ゲノム上の位置に偏った発現制御はないように見えます。siRNAによって抑制したERR alphaの結合モチーフ配列はWikipediaで見つけることができたので、これを転写開始位置近傍に持つ遺伝子を抽出することができました。これと、3種すべてのsiRNAで発現抑制された遺伝子リストとの重複を調べると、一つの遺伝子に辿り着きました。しかし、上述のとおり、３種すべてで抑制されている必要はないかもしれないですし、発現差解析の閾値がきつすぎだった可能性もあります。これらの条件を緩めることで、ERR alphaが直接制御する遺伝子の候補はもう少し広がるでしょう。
9 – おわりに
このチュートリアルでは、遺伝子発現量のデータを解析する手順を一通り紹介していますが、このやり方がどんな場合でも有効というわけではありません。実際のデータ解析では、データの特徴と研究目的に合わせた適切な判断が随所で必要になります。知識および経験不足のため適切な手法かどうかを判断することが難しい場合は、どうぞお気軽にデータ解析サービスについてお問い合わせください。