RNA-Seq解析では、FASTQファイルの品質確認、アライメント、Gene Countsの取得、正規化、発現変動遺伝子（DEG）の検出という一連の手順が広く使われています。Galaxy を使えば、RNA-Seq解析でよく利用されるツールをGUI上で順番に実行できます。コマンドライン環境を自分で構築しなくても、標準的な解析手順をたどれるのが大きな利点です。
では、こうした標準的な手順に従って各ツールを正しく実行すれば、最終的な解析結果が信頼できるかどうかまで確認できるのでしょうか。今回はGSE173789のデータを使って検証します。このデータセットについては、以前のケーススタディで本来の解析目的とは異なる要因によって生じたと強く疑われる発現変動が含まれていることを確認しています。まずはこのケーススタディ、または下記のショート解析1の動画で、このデータセットがどのような特徴を持っているのかを確認してから、この先を読み進めることをお勧めします。
このような明確な歪みを含んだデータを、Galaxyの標準的な手順に従って解析したとき、各工程で何が確認でき、どこまで結果の妥当性を判断できるのかを見てみました。解析手順については、Galaxy Training Network（GTN）のRNA-Seq reads to countsおよびRNA-seq counts to genesを参考にしています。
Galaxyでは、まず「解析を最後まで通すこと」に注意が向きやすい
Galaxyを実際に使って強く感じたのは、RNA-Seq解析が一つの統合ツールの中で完結するのではなく、複数のツールを順番につないでいく作業だということです。FASTQをFastQCへ渡し、その結果をMultiQCでまとめ、paired-endデータをHISAT2へ渡し、得られたBAMをfeatureCountsへ渡し、さらにGene Countsをlimmaで解析できる形式へまとめていきます。
実際には、それぞれのステップでかなり細かい問題が起こります。今回も、解析に入る前のFASTQ取得から問題がありました。Galaxy上でSRA accessionからFASTQを取得しようとしたとき、fasterq-dumpではエラーが発生し、利用者側からは原因を簡単に特定できませんでした。データ側の問題なのか、NCBI側なのか、Galaxy側なのか、ネットワークや一時的なサービス障害なのかを切り分けるのは容易ではありませんでした。この問題は、fastq-dumpに切り替えることで回避できました。
その後も、FastQCの結果をMultiQCへ渡すためにcollectionの形を整えるなど、ツール間でデータを受け渡すために形式を変えたり、整えたりする細かな作業がありました。HISAT2の実行では、2つのサンプルがsignal 9で終了する問題を回避する方法が分からず、この2サンプルを解析から除外することにしました。今回の目的ではそのような対応も可能でしたが、論文等の発表が目的の場合は、この問題を何とかして解決する必要があったでしょう。
こうした問題は、解析を先へ進めるためには一つずつ解決しなければなりません。そのため実際の作業では、どうしても注意が「なぜこのジョブが失敗したのか」「どの形式なら次のツールへ渡せるのか」「どうすれば最後まで処理を通せるのか」という方向へ向きます。エラーを解決して次のステップへ進み、最終的にDEG解析まで到達すると、それだけで一連の解析を無事に終えたような感覚も生まれます。
現在ではAIの助けが得られるので、このような問題への対応は以前とは比較にならないほど簡単になっています。それでも、解析を最後まで正常に実行することに注意が向きがちで、得られた結果が信頼できるかどうかを判断する難しさ自体は変わりませんでした。
FastQCとMultiQCには多くの情報があるが、除外判断まではできない
最初にすべてのFASTQについてFastQCを実行し、MultiQCで結果をまとめました。Sequence Counts、Sequence Quality、Duplication、Adapter Contentなど、多くのQC情報を一覧で確認できます。今回のデータでも、サンプル間でread数にはかなり大きな違いがあり、一部のサンプルではDuplicationやAdapter Contentにも他のサンプルとは明らかに異なる傾向が見られました。
もともと、このデータにはそのようなばらつきがあり、それと発現プロファイルが強く関連しているという予備知識があったので注意して見ていました。しかし、何の事前知識もない状態で同じレポートを見た場合、それだけを根拠にサンプルを除外すべきかどうかを判断することはできません。RNA-SeqではDuplicationが高いこと自体が必ずしも異常を意味するわけではなく、read数にもある程度のばらつきがあります。どの程度の違いから最終結果へ重大な影響を与えるのか、明確な境界が示されているわけではありません。
経験豊富な解析者であれば、QC情報から何らかの異常があるかもしれないと察知することは可能です。注意深い解析者であれば、気になるサンプルの番号をメモして、その後の工程でも追跡するかもしれません。しかし、初めて解析する人であれば、ほとんど気に留めることなく次の工程へ進んでいくことも十分に考えられます。
サンプル名がIDだけだと、レポートを直感的に理解しにくい
今回の解析でもう一つ感じたのは、サンプル名の付け方が結果やレポートの理解に大きく影響するということです。Galaxy内でSRA accessionからFASTQを取得すると、今回の場合はSRR accessionがそのままdataset名として使われました。そのため、MultiQCやfeatureCounts、後のMDSなどのレポートを見ても、SRR14409217やSRR14409222といった番号だけでは、それがMultipleSclerosisなのかHealthyControlなのかを直感的に判断できません。
もちろん対応表を見れば確認できます。しかし、解析工程をまたいで複数のレポートを比較するたびに対応表へ戻る必要があり、それだけでも判断の負担は大きくなります。特に今回のように、read数、Count分布、assignment率、MDS上の位置、DiseaseStateの関係を確認したい場合には、サンプル名に意味が見えないことが解析結果の理解を難しくします。
可能であれば、解析の初期段階で、元のIDを残しつつ主要な実験条件をsample labelへ含めておく方がよいでしょう。例えば、
SRR14409217_MS_sample01
のようにしておけば、各レポートの理解しやすさは格段に向上するでしょう。
ただし、今回のようにGalaxyに組み込まれたfastq-dumpなどを使って多数のSRA Runをまとめて取得した場合、生成されたdataset名を後から一つずつ変更するのはかなり煩雑です。それでも、解析前にサンプルIDと実験条件の対応を整理し、可能であれば早い段階で分かりやすいsample labelを設定しておくことが重要だと感じました。
featureCountsでも異常の手掛かりは見える
featureCountsを実行した後、MultiQCでassignment率を確認しました。多くのサンプルではおよそ55～64%程度でしたが、一部ではかなり低い値になっていました。しかし、やはりここでもfeatureCountsのassignment率だけで解析から除外すべきかどうかを判断することはできませんでした。assignment率が低いからといって、必ずしも後の発現量分布で同じ問題を示すわけでもないからです。
ここでも、レポートには注意すべき情報が現れてはいます。しかし、この段階で特定のサンプルを除外する判断をするのは難しく、そのままGene Countsを作成して次の統計解析へ進むのが自然でしょう。
そのままlimma-voomまで到達してしまった
featureCountsで得られたGene Countsを結合し、limma-voomによる発現変動解析を行いました。比較したのはMultipleSclerosisとHealthyControlです。low-count genesについてはCPMに基づくフィルタリングを行い、TMMによる正規化とvoomによるmean-variance modelingを行いました。
limmaのレポートにはMDS plot、mean-variance trend、model fit、MD plot、Volcano plotなどが表示されます。これらを見る限り、解析はきちんと実行されています。mean-variance trendも自然に見え、Volcano plotには多数の発現変動遺伝子候補が表示されます。adjusted p-valueを基準にDEGを抽出すれば、そのままGO解析やPathway解析へ進むこともできます。
MDS plotを見ると、MultipleSclerosisとHealthyControlも一定程度分離しているように見えます。DiseaseStateだけを色分けしてこの図を見れば、「MultipleSclerosisとHealthyControlでは遺伝子発現プロファイルが異なっている」と解釈するのは自然です。標準的な解析を最後まで実行し、診断プロットにも大きな違和感がなく、さらにDEGまで得られているため、ここまで来ると解析結果を信頼したくなります。
しかし、解析結果を元のGene Countsと合わせて見ると、別の構造が見える
そこで、featureCountsで得られたGene Countsと解析結果をSubio Platformでまとめて可視化してみました。元のGene Countsの分布をサンプルごとに見ると、一部のサンプルでは発現量全体が他のサンプルよりかなり低いことが分かります。さらにSampleGroupで色分けすると、doubtやcloseToDoubtに分類されているサンプルの多くが、この低い分布側に位置しています。これらは、上述のケーススタディで注意を喚起したサンプル群です。
正規化後の値を見ると、このサンプル間の分布差はかなり小さくなります。もちろん、RNA-Seqではlibrary sizeやRNA compositionが異なるため、Gene Countsをそのまま比較することはできず、正規化は必要です。しかし今回重要なのは、補正後に見えている差が生物学的な差なのか、それとも元々存在していた技術的な差と関係しているのかという点です。
FASTQのQCレポートで見られたread数のばらつき、featureCountsで見られたassignment率の違い、元のGene Countsの分布、正規化によってどのように値が変化したのか、そしてMDSやDEGの結果をまとめて見ることで、初めてそれぞれの情報の関係が見えてきます。少なくとも、最終的なDEGだけを見ていたのでは気づけなかった問題を、解析全体の中で検討できるようになります。
limmaが間違っているわけではない
私たちは、Galaxyやlimma-voomを批判したいわけではありません。limmaは、与えられた発現データと実験デザインに基づいて統計解析を行うツールです。入力データが解析対象として妥当であり、指定された比較を行ってよいという前提のもとで、library sizeやmean-variance relationshipなどを考慮しながら群間差を評価します。
しかし、limmaに到達したデータが、本当にその前提を十分に満たしているとは限りません。FastQC、HISAT2、featureCountsという標準的な工程を通過したからといって、最終的な群間比較に適したデータになっていることまで保証されるわけではないからです。
言い換えれば、limmaは与えられたデータを統計的に解析しますが、そのデータをそもそも比較してよいのかまでは判断しません。その判断は解析者に残されています。
各ツールのレポートだけでは、解析全体は見えない
FastQCはFASTQの品質を評価し、featureCountsはreadの遺伝子への割り当てを評価し、limmaは与えられた発現行列とdesign matrixに基づいて統計解析を行います。どのツールも、それぞれの目的の範囲では必要な情報を返しています。
しかし、FastQCでread数の少なかったサンプルがfeatureCountsでどの程度のCountsを持ち、limmaのMDSでどこに位置し、さらにどのDiseaseStateに属しているのかまでをまとめて示してくれるわけではありません。今回も、一つ一つのレポートでは解析方針を変えるほど強い異常には見えなかった情報が、同じサンプルについてまとめて見ることで、一つの傾向として見えてきました。
そのため、重要なのはQC情報を増やすことではなく、read数、assignment率、Gene Countsの分布、正規化後の変化、MDS、DEG、さらに実験条件までを、同じサンプルについてまとめて確認できることです。サンプル名がSRR番号だけでは理解しにくかったのも、このためです。
それでも、解析者は判断しなければならない
すべての情報をまとめて見ても、サンプルを解析に含めるべきか、除外すべきかを一意に決められるとは限りません。RNA-Seqデータにはもともと大きなばらつきがあり、read数が少ない、Duplicationが高い、assignment率が低い、MDSで離れているといった特徴について、どこからを異常とするかという普遍的な境界はありません。
特に、技術的な違いとDiseaseStateが重なっている場合には、観測された差が生物学的なものなのか、技術的なものなのか、あるいはその両方なのかを簡単に判断することはできません。それでも、最終的にそのデータを解析に含めるのか、結果をどこまで信頼するのかを解析者は判断しなければならないのです。
だからこそ、その判断に必要な情報は、別々のレポートとして断片的に提示されるべきではありません。read数、assignment率、Gene Countsの分布、正規化による変化、MDSやDEG、さらにサンプルIDと実験条件など、まとめて確認するための手段がどうしても必要なのです。
パイプラインを動かすことと、解析することは別
以前、RNA-Seq解析で、ヒートマップとDEG解析に異なる値を使ってもよいのか？では、異なるツールが異なる値を使って解析しているにもかかわらず、その違いが十分に意識されないまま、それぞれの結果が同じ現象を別の方法で示したものとしてまとめられることの問題を取り上げました。実際には、それぞれのツールが扱っている情報は本質的に異なります。
今回Galaxyで見えてきた問題はそれとは別のものです。FastQC、featureCounts、limmaなど、それぞれの工程には重要な情報が含まれていても、一つ一つを単独で見れば、サンプルを除外したり解析方針を変更したりするほど強い材料には見えないことがあります。そのため、断片的な情報として見過ごされ、そのまま次の工程へ進んでしまいます。しかし、read数、assignment率、Count分布、正規化前後の変化、MDS、DEG、さらにDiseaseStateなどの実験条件まで、同じサンプルについてまとめて見ると、それまで弱かった情報が一つの傾向として見えてきます。
これらは別の問題ですが、原因は共通しています。つまり、RNA-Seq解析では、それぞれ単独では完成されたツールをつなげて一つの解析を組み立てていることにあります。個々のツールの結果を順番に確認していくだけでは、解析全体を通して見たときに何が起きているのかを捉えることはできません。問題は個々のツールではなく、複数の工程を通して得られた情報を、解析全体として統合して解釈するという視点が、現在のパイプライン型の解析では十分に意識されていないことにあります。
この問題はGalaxyに限ったものではありません。BioconductorやPythonなどで複数のツールやパッケージを組み合わせて解析する場合にも同じです。これまでRNA-Seq解析では、決められたパイプラインを正しく実行し、各ツールの設定やエラーに対処すること自体が、解析者に求められる重要なスキルでした。しかし、それだけでは解析とは言えません。
信頼性の高いツールで構成されたパイプラインを最後まで動かせたからといって、得られた解析結果が信頼できるとは限りません。複数の工程を通して得られた情報を、解析全体として統合して解釈すること。それこそが解析者の役割です。
関連リソース： RNA-Seq解析の一連の流れを実際のデータを使って学びたい方は、RNA-Seqデータ解析チュートリアルをご覧ください。解析の各段階で得られた結果を可視化し、元のデータと照らし合わせながら確認・検証するには、Subio Platformを利用できます。
