RNA-Seq解析では、発現差解析で抽出した遺伝子を、 ヒートマップとして可視化することがよくあります。 その際、元の発現値ではなく、 遺伝子ごとにZ-scoreへ変換した値や、 変動を0～1の範囲へscalingした値が 使われることがあります。
Z-score化やscalingは、 ヒートマップを描く前処理として広く利用されているため、 RNA-Seqのヒートマップを作成するために 必要な処理だと思われることがあります。 しかし、 実際には必須の処理ではありません。
この記事では、Z-scoreやscalingによって 何が保持され、何が失われるのかを整理します。 そのうえで、一般的なRNA-Seqデータセットでは、 row scalingが必要ない場合が多いだけでなく、 ほかの解析結果との整合性を損なう可能性があることを説明します。 さらに、row scalingが必要なケースについても検討します。
用語の整理
ヒートマップの各行について、 平均値を引くなどして基準位置をそろえる処理を centeringといいます。 これにより、遺伝子ごとの発現量の違いを取り除き、 平均からの増減や遺伝子間の連動性に注目できます。 centeringしても、遺伝子ごとの変動幅の違いは保持されます。 RNA-Seqで対数変換した値をcenteringする場合、 発現量の差分ではなく、 各遺伝子の基準値からの相対的な変化（fold change）に 注目していることになります。
これに対して、 標準偏差や値の範囲などを用いて、 各行の変動幅までそろえる処理を scalingといいます。 ヒートマップの各行に対して個別にscalingを行うことを、 この記事ではrow scalingと呼びます。
row scalingの方法には、Z-scoreへの変換や、 最小値と最大値を用いるmin–max scalingなどがあります。
Z-scoreへの変換では、 各遺伝子の平均値を引いてcenteringした後、 その遺伝子の標準偏差で割ります。 その結果、各遺伝子の平均は0、 標準偏差は1になります。 したがって、もともとの変動幅が小さい遺伝子も大きい遺伝子も、 ヒートマップ上では同程度の変動を持つように表示されます。
min–max scalingでは、 各遺伝子の最小値を引いた後、 最大値と最小値の差で割ります。 その結果、各遺伝子の最小値は0、 最大値は1になります。 例えば、100から200までの変動も、 100から10,000までの変動も、 0から1の範囲へ変換されます。 この記事では、この変換を 0～1 scalingと表記します。
Note: VSTやrlogをheatmapに使用する場合
なお、RNA-Seq解析でよく使われるVSTやrlogは、 ここでいうrow scalingとは異なります。 VSTやrlogは、Countsの平均値と分散の関係を緩和し、 低Counts領域の大きなばらつきが 解析結果に過度に影響することを抑えるための変換です。 これらを適用してからclusteringやheatmapの描画を行うと、 低Counts領域の遺伝子のばらつきが 相対的に目立ちにくくなります。 そのため、低Counts領域を十分に除外したデータでは、 VSTやrlogと単純なlog変換との差は小さくなります。
VSTやrlogをheatmapに使用することが妥当かどうかは、 DEG解析に用いた方法と値によっても異なります。 DESeq2では、Countsを用いて発現差を検定し、 VSTやrlogは主にclusteringや可視化のために用いられるため、 DEGとして抽出された遺伝子のサンプル間パターンを確認する値としては、 実用的な近似になります。 一方、edgeRの結果を確認する場合には、 TMM正規化後のlogCPMなど、 edgeRの正規化に対応した近似値を表示に用いる方が自然です。 また、log変換した正規化値やTPM、FPKMを用いて t-testやノンパラメトリック検定を行った場合には、 heatmapにも同じ値を用いる方が、 検定結果との整合性は高くなります。
ただし、このような可視化用の近似値を用いると、 元のGene Countsや サンプルごとのダイナミックレンジとの関係は 見えにくくなります。 したがって、VSTやrlogが 常に適切な可視化用の値であるとは限りません。 DEG解析に用いた方法や値との関係だけでなく、 元データのどの特徴を確認したいのかも考えたうえで、 表示に用いる値を選ぶ必要があります。
以上をまとめると、 VSTやrlogは、 DEGとして抽出された遺伝子のサンプル間パターンを 分かりやすく提示する目的には適しています。 一方で、元のGene Countsとの関係が見えにくくなるため、 解析結果がどのように生じたのかを検証する目的には、 適さない場合があります。
なぜ遺伝子ごとに値を変換するのか
遺伝子ごとに発現量は大きく異なります。 たとえば、ある遺伝子がgene counts 1～10の範囲で変化し、 別の遺伝子が 100～1000の範囲で変化しているとします。 すべての遺伝子に共通するカラースケールを使うと、 色の違いは発現量の差を表すことになり、 個々の遺伝子の発現変動が見えにくくなります。
この問題を避けるために、centeringは必要です。 しかし、row scalingを行って、 遺伝子ごとの変動幅を均一化することまでも必要でしょうか？
振幅は生物学的解釈に必要な情報である
RNA-Seqデータにおける変動幅は、 単なる数値上の違いではありません。
遺伝子発現の変化量は、 応答の強さや生物学的な影響の大きさを考えるための 重要な情報です。 たとえば、p値が同じであっても、 1.2倍しか変化していない遺伝子と、 10倍変化している遺伝子とでは、 表している生物学的現象が全く異なることがありえます。
振幅は他の解析結果の妥当性を検討するために必要な情報である
また、振幅の情報は、 発現差解析やほかの解析結果の妥当性を確認するためにも必要です。 たとえば、DEGとして抽出された遺伝子をヒートマップに表示する場合、 少なくとも次の点を確認する必要があります。
- 平均値で見える群間差が、個別のサンプルでも確認できるか
- 反復サンプルで変化方向が再現されているか
- 一部のサンプルだけが差を作っていないか
- 同一実験区内で、発現変動の幅が異なるサブグループはないか
しかし、row scalingを行えば変化量の違いが失われ、 その結果、上記のような違いが見えなくなる可能性があります。
さらに、実験区が複数ある場合に、 各実験区の繰り返しサンプルの平均値をヒートマップにしてしまうと、 各群の群内のばらつきが見えなくなります。 したがって、 平均化とrow scalingを組み合わせたヒートマップは、 DEG解析の妥当性を確認するために必要な情報を 二重に減らしてしまうことになります。
振幅を均一化する必要がある場面は限定的である
それでは、振幅を均一化することが 合理的になるのはどのような場合でしょうか。
一般的な統計解析では、 身長と体重のように、 単位も数値範囲も異なる特徴量を 同じ多変量解析へ入力するような場合、 数値範囲の大きい変数だけが結果を支配しないように、 各特徴量の尺度をそろえる必要があります。 これが、scalingが必要とされる最も一般的な理由です。
しかし、一つのRNA-Seqデータセット内では、 各遺伝子の発現値は、 基本的に同じ測定法、同じ単位、同じ解析手順によって得られています。 遺伝子ごとの変動幅の違いは、 発現量や応答強度、測定の安定性などに関する情報です。 そのため、このような重要な情報を捨ててまで row scalingを行うメリットは、通常ありません。
ただし、RNA-Seqにおいてscalingが必要になる可能性があるのは、 たとえば次のような場合です。
- 同じRNA-Seqでも、独立した複数のRNA-Seqデータセットを統合する場合
- 同じトランスクリプトームでも、RNA-SeqとマイクロアレイやqPCRデータなど測定手法の異なるデータセットを統合する場合
- 異なる種類のオミクス（例えばプロテオームやDNAメチル化など）のデータセットを統合する場合
- 発現データと臨床検査値などを同じ解析へ入力する場合
ただし、この場合でも、 単純にZ-scoreへ変換すれば 測定法間の違いやバッチ効果が解決するわけではありません。 まず各データについて適切な正規化や品質確認を行い、 その後の解析目的とアルゴリズムに応じて、 必要であれば尺度をそろえるという順序が必要です。
また、機械学習の前処理としてscalingを使う場合も、 後段のアルゴリズムが特徴量の尺度に影響されるかどうかによって 必要な処理方法は異なります。
したがって、 振幅を均一化することは一般的な前処理ではなく、 解析目的と後段の手法に依存する、条件付きの選択肢です。
GSE95304のFig. 5Fを再解析する
ここから、公開RNA-SeqデータGSE95304を用いた研究の Fig. 5Fを見ていきます。 このデータは、別の記事でも例として使っていますが、 この論文の解析の進め方やロジック構成を批判するのではなく、 一般的なRNA-Seq解析にも共通する問題を考えるための 具体例として取り上げます。
この研究では、 Fra-1、MYC、E2F1、TP53などの転写因子を shRNAによってknockdownし、 遺伝子発現への影響を調べています。 論文では、FPKMを用いた発現解析から 295遺伝子を抽出し、 Fig. 5Fにヒートマップとして示しています。
公開されたFPKMデータに、 Methodsに記載されたフィルターと検定を適用したところ、 23,615遺伝子から12,239遺伝子へ絞り込まれ、 FDR < 0.1を満たす295遺伝子を再現できました。 このうち、157遺伝子はknockdown群で低下し、 138遺伝子は上昇していました。 したがって、 Fig. 5Fに使われた遺伝子群の選択過程を、 公開データからおおむね再現できたと考えられます。
295遺伝子は、各shRNAで有意だった遺伝子ではない
この295遺伝子の意味を理解するには、 どのような比較によって抽出されたのかを確認する必要があります。 295遺伝子は、 各shRNAについて個別に発現差解析を行い、 その和集合を取ったものではありません。 8個のshSCRサンプルをcontrol群とし、 Fra-1、MYC、E2F1、TP53のknockdown各2サンプル、 合計8サンプルをknockdown群として、 8対8の二群間比較が行われています。 したがって、この検定によって示されるのは、 4種類のknockdownをまとめた群と、 shSCR群との間に平均的な発現差があった ということです。
Fig. 5Fは、検定とは別の問いに答えているように見える
Fig. 5Fでは、 295遺伝子がshRNA条件ごとに分けて表示されています。 しかも、各条件の反復サンプルは平均化され、 さらに各遺伝子の値には0～1 scalingが適用されています。 ただし、平均化とscalingが どの順序で行われたのかは明確に記載されていません。 いずれの順序で処理されたとしても、 Fig. 5Fからは、 各shRNAによって個別に影響を受けた遺伝子を集めたところ、 それらの多くが複数のshRNA条件で 連動して変化していたかのような印象を受けます。 発現差解析が答えた問いと、 Fig. 5Fが読者に提示しているように見える内容は、 完全には一致していません。
この発現差解析によって捉えられた群間差を ヒートマップで確認するのであれば、 発現差解析と同じ値、 または可視化に必要な最小限の対数変換とcenteringだけを加えた値を用いて、 発現差解析に使われたものと同じサンプルを 個別に可視化する必要があります。
この図では、発現差解析に使われたFPKMを対数変換し、 遺伝子ごとにcenteringした値を用いています。 発現差解析に使われた16サンプルを、 検定と同じく8個のcontrol群と8個のknockdown群として配置し、 平均せずに個別に表示しています。 上段の「DEG group」は、 発現差解析で用いられた二群の区分を示します。 その下には、各サンプルのshRNA条件と反復番号を示しています。 また、遺伝子ごとの変動幅を保持するため、 row scalingは行っていません。
shFra-1、shp53、shMycでは、 2つの反復サンプルが比較的似た発現パターンを示しています。 一方、shE2F1では、batch 2のサンプルで発現パターンが異なっています。 shE2F1でも、ほかの3種類のshRNAと同様の変化が生じている可能性はありますが、 今回の2サンプルだけからは判断を保留する必要があります。
また、controlであるshSCR群でも、 すべてのサンプルが均一な発現パターンを示しているわけではありません。 特に、batch 1とbatch 2の間で、 発現パターンに違いが見られます。
このように個別サンプルを表示すると、 抽出された295遺伝子の群間差には、 shSCRのbatch 2の4サンプルと、 shE2F1のbatch 2サンプルを除く7個のknockdownサンプルとの違いが 強く反映されていることが分かります。
そして重要なことは、 Fig. 5Fの代わりにこの図を用いたとしても、 複数のshRNA処理によって共通して変化する遺伝子が存在するという 論文の中心的な考え方を支持することに変わりはないということです。 むしろ、個別サンプルのばらつきを示したうえで、 どのknockdown条件で共通性が確認でき、 どの条件では判断を保留すべきかを区別できる情報を 提示する図になります。
解析結果を関連付けるためには、値の一貫性が重要である
ここまで見てきた問題は、 特定の可視化方法だけに限られるものではありません。 より一般的には、 複数の解析結果をどのような根拠で関連付けて解釈するのか という問題です。
RNA-Seq解析では、 PCA、クラスタリング、ヒートマップ、 発現差解析、エンリッチメント解析など、 複数の結果が一つの流れとして示されます。 しかし、それぞれの結果が 異なる値、異なるサンプル構成、 異なる遺伝子の選択条件から作られていれば、 同じデータセットを解析した結果であっても、 互いに同じ内容を示しているとは限りません。
複数の解析結果を関連付けて解釈するためには、 可能な範囲で、 同じ値、同じサンプル、同じ遺伝子群を使う方が明確です。 解析間で共通する情報が多いほど、 一つの結果で見られた特徴が 別の結果にも表れているかを直接確認できます。
一方、解析目的や手法上の理由から、 異なる値や異なるデータ処理を使う必要がある場合もあります。 その場合には、 それぞれの処理によって何が保持され、 何が変化し、 何が失われたのかを明確に認識する必要があります。 また、結果を提示する際にも、 どの値に、どのような処理を加えて作成した結果なのかを 明示すべきです。
このような確認は、 Rのツールを個別に組み合わせて作成した解析パイプラインや、 複数の解析ツールを自動的に実行するパッケージソフトでは、 必ずしも容易ではありません。 各解析がどの値から作られ、 途中で利用者が意図していない変換や、 把握していない変換が行われていても、 それに気づきにくいためです。
そのため、一連の解析結果が自動的に作成される場合ほど、 個々の解析が正しく実行されたかだけでなく、 各結果がどの値から作られ、 結果どうしをどこまで関連付けて解釈できるのかに 注意する必要があります。 この問題については、 ヒートマップで可視化された値とDEG解析に使われた値は、 異なるかもしれない でも詳しく説明しています。
