RNA-Seq解析の結果を、どこまで検証できるのか
以前の記事「403：RNA-Seqデータの再解析」では、多発性硬化症（Multiple Sclerosis：MS）のRNA-Seqデータを再解析し、サンプル間の発現値の分布に無視できない違いがあることを確認しました。一部のMSサンプルでは、多数の遺伝子の発現値が高くなっており、それらが発現上昇遺伝子（Up DEG）の判定に影響していました。また、Healthy Controlのサンプルも、発現パターンの異なる二つの集団に分かれる傾向が見られました。このようなデータでは、DEG解析の結果をそのまま受け入れるのではなく、どのようなサンプル間の違いによって結果が生じているのかを調べることが重要です。
これまで、Galaxy やRaNA-Seq を使って、同じデータを解析し、こうした問題をどこまで発見・検証できるかを調べてきました。今回は、RNAseqChefを試してみます。
RNAseqChefは多様な解析機能を備えている
RNAseqChefは、豊富な可視化および解析機能を備えていると謳っているので、解析結果をさまざまな角度から調べ、データに問題が見つかった場合にも、その原因を詳しく検証できるのではないかと期待しました。しかし、実際に使ってみると、基本的な操作方法や結果の提示方法はRaNA-Seqとよく似ていました。データを入力すると一連の解析が自動的に実行され、結果がレポートとして表示されます。
たしかに豊富な解析機能があるのですが、複数のグラフを行き来しながら、特定のサンプルや遺伝子について問題の原因を掘り下げることのできるユーザーインターフェースは備えていませんでした。解析・可視化機能が豊富であることと、解析結果を自由に検証できることは、同じではないようです。
Gene Countsを入力すると、自動的に解析結果が得られる
RNAseqChefでは、解析対象データセットのGene Countsとサンプル情報を記載したテキストファイルをアップロードします。ただし、遺伝子IDにはEnsemblかGene Symbolのみ受け付け、NCBI Gene IDは対応していません。GalaxyやRaNA-Seqと異なり、FASTQファイルの処理は含まれていません。今回は、できるだけ同じ条件で比較したかったので、Galaxyで取得したGene Countsをアップロードしました。
DEG解析は、手法としてDESeq2、edgeR、EBSeqを選択でき、FDRの補正方法、Fold Change、FDR、BaseMeanなどの閾値も設定できます。今回は、すべてデフォルトの設定で解析を実行しました。解析を実行すると、PCA、クラスタリング、MA plot、ヒートマップ、DEG解析の結果などが自動的に生成されます。
前のGalaxyやRaNA-Seqでの検証との大きな違いは、解析を始めるための入力データです。RNAseqChefではGene Countsを入力するため、総リード数や遺伝子にアサインされたリード数など、FASTQからGene Countsを生成する過程のQC情報は得られません。これは入力データの違いによるものです。
一方、入力されたGene Countsだけでも、サンプルごとの発現値の分布を比較することは可能なはずです。しかし、今回使用したPair-wise DEGモードには、そのための図が用意されていませんでした。今回のデータで問題となった発現分布の歪みを調べるための、重要な手がかりが得られないのです。
PCAやヒートマップには、問題のあるデータ構造が現れている
RNAseqChefのResult overviewには、PCA、MDS、階層的クラスタリング、MA plot、ヒートマップなどが表示されます。今回のデータでも、403の記事で指摘した特徴的な構造が確認できました。
PCAでは、Healthy ControlとMSのサンプルは明瞭に分離せず、一部のサンプルが大きく外れています。また、ヒートマップには、一部のMSサンプルで多数の遺伝子の発現値が高くなっていることに加え、Healthy Controlのサンプルも異なる発現パターンを示す集団に分かれていることが現れています。この二段構えの構造は、原著論文のFig. 4にも現れています。つまり、RNAseqChefの解析結果にも、データを慎重に検討する必要があることを示す手がかりは含まれています。
PCAで外れたサンプルをどう解釈するか
PCAで他のサンプルから大きく外れたものが見つかっても、それだけを理由に解析から除外するという判断はできません。PCAの結果からだけでは、その違いが生物学的な違いを反映したものではないと断定することができないからです。
特に今回のデータでは、MSサンプルの一部に特徴的な発現パターンが見られます。これを見て、MSには異なるサブタイプが存在するのではないか、という仮説に進むことも考えられます。MSは生物学的に不均一な疾患であり、サブタイプの存在を考えること自体は不自然ではありません。
しかし、403の記事で確認したように、今回のサンプル間の違いには、Gene Countsの分布の歪みが関係しています。そこまで調べることができれば、観察された構造の少なくとも一部がGene Countsの分布の歪みに由来することを確認でき、生物学的なサブタイプとして解釈することが適切なのか、慎重に判断するための根拠が得られます。しかし、このような調査に必要な図やユーザーインターフェースをRNAseqChefは提供していません。
また、Healthy Controlのサンプルも二つの集団に分かれるように見えることも、解析結果の妥当性を検討するうえで無視できません。これについては解析者がヒートマップを見てどう判断するかに委ねられますが、これらのサンプルを深追いする機能はありません。
RNAseqChefでは、表示されたグラフは基本的に独立しており、PCAで気になったサンプルを選択して、別のグラフで対応するデータを確認することはできません。例えば、大きく外れたサンプルが見つかっても、そのサンプルでは発現値が全体としてどのように分布しているのかを、そのまま追跡することはできませんでした。
Z-scoreヒートマップでは、遺伝子間の変動幅を比較できない
RNAseqChefのヒートマップは、Z-scoreで表示されています。遺伝子ごとに平均値を引くCenteringは、遺伝子間の発現量の違いを取り除き、サンプル間の発現変動を比較するために有用です。しかし、Z-scoreではさらに標準偏差で割ることで、遺伝子ごとの変動幅がそろえられます。（詳しくは、「RNA-Seqのヒートマップで、Z-scoreを安易に使うべきでない理由」をご覧ください。）
今回のように、一部のサンプルにおける発現値の分布の歪みが大きな変動を生み出している場合、その結果の解釈に大きな影響を与えることが予想されるため、特に注意が必要です。しかし、RNAseqChefではこのZ-scoreによる表示が固定で、変更できません。
GOI profilingでは、個々の遺伝子の発現分布を確認できる
RNAseqChefには、興味のある遺伝子を指定して解析結果を確認するGOI profilingという機能があります。遺伝子名を入力すると、Volcano plot上での位置やヒートマップ、正規化Countsを使ったBox plotなどが表示されます。これは便利な機能です。MA plotやVolcano plotでは個々のサンプルの発現値のばらつきが分からないため、Box plotでサンプルごとの値を確認できる点は評価できます。
しかし、Box plot上で大きく外れた値が見つかっても、その点がどのサンプルに対応しているのかを確認することはできません。さらに、そのサンプルでは他の遺伝子も同様に高い値を示しているのか、発現値の分布全体が他のサンプルと異なっているのかも調べることができませんでした。
DEG解析からエンリッチメント解析まで自動的に進む
RNAseqChefでは、DEG解析で抽出された遺伝子が自動的にエンリッチメント解析に使用されます。この点もRaNA-Seqと共通しています。今回のように、DEG判定にサンプル間の発現分布の歪みが影響している可能性がある場合、DEGの妥当性を十分に検証できないまま、次の解析へ進んでしまうことになります。エンリッチメント解析を自動実行できること自体は便利ですが、入力する遺伝子リストが適切かどうかは、解析者自身が判断しなければなりません。
一方、RNAseqChefのエンリッチメント結果では、グラフの下に表示されるテーブルを右にスクロールすると、各Gene setに含まれる該当遺伝子の名前を確認できます。そのため、異なるGene set間でどの遺伝子が重複しているかを、Webブラウザのページ内検索機能でも調べられます。グラフを見るだけで終わらず、結果に寄与した遺伝子まで確認できる点は評価できます。
RaNA-SeqとRNAseqChefを使い分けるなら
実際にRNAseqChefを使ってみると、基本的な操作方法や結果の提示方法はRaNA-Seqとよく似ていました。どちらも、データを入力すると標準的な解析が自動的に進み、PCA、クラスタリング、DEG解析、エンリッチメント解析などの結果が提示されます。もちろん機能を細かく比較すれば様々な違いはあります。
しかし、ユースケースで考えた場合、今回検証したような自動解析を利用するのであれば、FASTQを持っている場合はRaNA-Seq、Gene Countsを持っている場合はRNAseqChefを使うことになります。さらに、3群以上の多群間比較を行うのであれば、RaNA-SeqからGene Countsをダウンロードして、RNAseqChefに持っていくこともできるでしょう。
解析結果が得られることと、その信頼性を検証できることは違う
今回の検証では、RNAseqChefを使ってDEG解析からエンリッチメント解析まで簡単に実行できました。PCAやヒートマップには、403の記事で指摘した問題のあるデータ構造も現れています。また、GOI profilingを使えば、興味のある遺伝子の発現分布まで確認できます。しかし、そこで見つけた疑問をさらに掘り下げることはできませんでした。
解析結果が得られることと、その結果がどのようなデータの特徴から生じたのかを理解し、信頼性を検証できることは、別の問題です。
これはRNAseqChefに限った問題ではありません。Galaxy、RaNA-Seq、RNAseqChefでは解析の自由度や操作方法に違いがありますが、今回の検証では、解析結果に現れた問題の原因を調べ、その妥当性を判断することが、いずれのツールでも容易ではありませんでした。こうした問題の背景には、データ解析を学ぶ際に、「ツールを使って結果を出すこと」に比べて、「得られた結果が信頼できるかを検証すること」が十分に重視されてこなかった可能性も考えられます。
しかし、「Designing RNA sequencing experiments: A practical guide to reproducible gene expression analysis」では、RNA-Seq解析における再現性やデータ解釈の問題を指摘し、実験設計、品質管理、解析手法の選択などを通して、再現性があり、生物学的に意味のある結果を得ることの重要性を論じています。
重要なのは、PCAやクラスタリングなどの図を作ることではなく、そこに何が現れているのかを調べ、得られた結果が信頼できるものなのかを判断することです。
解析ツールが自動的に多くの結果を生成できるようになった今だからこそ、解析者自身がその結果を検証できることの重要性は、むしろ高まっているのではないでしょうか。
