この記事は、前回のケーススタディを前提とした続編です。
本記事では、前回観察されたサンプル構造について、
正規化方法や補正方法を変えた場合にどう見えるかを追加検討しています。
サンプル除外や補正に関する判断は、本記事だけで独立に行ったものではなく、
前回のケーススタディ
での可視化結果を前提にしています。
まず前回の記事を読んだうえで、本記事をご覧ください。
前回の記事では、RNA-SeqデータをCPM正規化して可視化したときに、 一部のサンプルが他のサンプルとは異なる発現プロファイルを示し、 クラスタリングでも特徴的なまとまりを形成していることを紹介しました。
このような結果を見ると、 「CPM正規化の影響なのではないか」 「別の正規化方法を使えば解決するのではないか」 と考えたくなります。
そこで今回は、同じデータに対して複数の正規化・変換・補正方法を試し、 前回観察されたサンプル構造がどの程度変化するのかを確認しました。 また、クラスタリング結果だけでなく、各サンプルのヒストグラム上で、 doubtサンプルに特徴的な遺伝子群がどのように見えているのかも確認しました。
ここで重要なのは、どの方法が唯一の正解かを決めることではありません。 正規化や補正を行った後に、データの見え方がどのように変わるのか、 または変わらないのかを確認し、 その構造が何を反映している可能性があるのかを考えることです。
正規化や変換を変えても、doubtサンプルの構造は残った
図中のカラーブロックでは、前回の解析で観察されたサンプル構造を示しています。 赤は、PCAやクラスタリングで他のサンプルとは明確に異なる位置に現れたdoubtサンプルです。 黄色は、PCA上ではdoubtサンプルと大部分のその他のサンプルの中間に位置し、 ヒストグラムで確認してもdoubtサンプル寄りの発現プロファイルを示したため、 close to doubtとラベル付けしたサンプルです。
緑は、doubtおよびclose to doubtのサンプルを除いた後、 残りのサンプルだけでクラスタリングを行ったときに検出されたクラスターです。 このクラスターには、Healthy controlサンプルの約半数が含まれていました。 青はそれ以外のサンプルで、残りの約半数のHealthy controlとDiseasedサンプルが含まれています。
TMM正規化後のCPMでも、doubtサンプルの構造は残った
まず、同じGene Countデータに対してTMM正規化を行い、 TMM正規化後のCPMデータを用いてクラスタリングを行いました。 TMM正規化は、RNA-Seq解析で広く使われている有用な正規化方法です。
しかし、このデータでは、TMM正規化後のCPMを用いても、 前回のCPM正規化データで見えていたサンプル構造は大きく変わりませんでした。 図1-aの最下部のカラーブロックは、前回のCPM正規化データに基づくクラスタリング結果を示しています。 赤で示したdoubtサンプルと、黄色で示したclose to doubtサンプルは、 TMM-CPM補正後も近い発現プロファイルを持つサンプル群としてまとまりました。
今回は、doubtサンプルで共通して高く見える遺伝子クラスターをhigh@doubt、 低く見える遺伝子クラスターをlow@doubtとして取り出しました。 図1-aでは、これらの遺伝子クラスターを灰色のオーバーレイで示しています。
DESeq2 VSTでも、doubtサンプルの構造は残った
次に、同じGene CountデータをDESeq2のvariance stabilizing transformation（VST）で変換し、 その値をSubio Platformに取り込んでクラスタリングを行いました。 VSTは、カウント値の平均と分散の関係を緩和し、 サンプル間の距離やクラスタリングを評価しやすくするための変換です。
VST値はすでに変換済みの値であるため、 Subio Platform側では追加のlog2変換を行わず、そのまま使用しました。
VST変換後のデータでも、 赤で示したdoubtサンプルと、黄色で示したclose to doubtサンプルは、 再び近い発現プロファイルを持つサンプル群としてまとまりました。 また、図2-aでは、図1-aで取り出したhigh@doubtおよびlow@doubtの遺伝子クラスターを 灰色でオーバーレイしています。 VST変換後でも、これらの遺伝子群はdoubtサンプルに特徴的な発現パターンとして確認できました。
Quantile normalizationでも、クラスタ構造は消えなかった
さらに参考比較として、Raw Gene Countデータに対してquantile normalizationを行いました。 Quantile normalizationは、サンプル間の分布を強制的にそろえる方法であり、 マイクロアレイ解析ではよく用いられてきました。
一方で、RNA-Seqデータでは、 生物学的に意味のある発現分布の違いまでそろえてしまう可能性があるため、 標準的な解決策としてではなく、 分布を強く補正した場合の参考として扱う必要があります。
今回は、Gene Countデータにquantile normalizationを行い、 その後log2変換とcentering（ratio to average）を行ってクラスタリングしました。 しかし、この強い補正を行っても、 前回のCPM正規化データで見られたクラスタ構造は大きく変わりませんでした。 図3-aでも、high@doubtおよびlow@doubtの遺伝子クラスターは、 doubtサンプルに関連する特徴的なパターンとして残っていました。
ヒストグラムで見ると、補正後もdoubtサンプルの特徴が残っている
Gene Countでは、doubtサンプルのヒストグラムが左に寄っている
クラスタリングだけでなく、各サンプルのヒストグラムを比較すると、 doubtサンプルの特徴がより分かりやすくなります。 図4では、Gene Count、TMM-CPM、VST、Quantile normalization後の分布を並べて比較しています。
Gene Countデータでは、doubtサンプルのヒストグラムが分かりやすく左側に寄っています。 これは、これらのサンプルで全体的にcount値のダイナミックレンジが狭くなっていることを示しています。
TMM-CPMでは、ヒストグラムの形を比較的保ったまま、 分布全体が横方向に移動しているように見えます。 このことから、TMM-CPMでは比較的線形に近い補正が行われていることが分かります。
一方、VSTでは分布の形そのものが大きく変わっています。 これは、VSTがカウント値の平均と分散の関係を緩和するための変換であり、 単純にヒストグラムを横にずらす補正ではないためです。
Quantile normalizationでは、基本的にサンプル間の分布形をそろえる方向に補正されます。 そのため、doubtサンプルのように、もともとの分布が他のサンプルと大きく異なり、 かつ該当するサンプル数が少ない場合には、 補正後の分布形も大きく変化します。
VSTおよびquantile normalization後のヒストグラムでも、 doubtサンプルでは左側が欠けたような形が残っていました。 このことは、もともとのダイナミックレンジの狭さが、 補正後の分布にも影響していることを示していると考えられます。
high@doubt遺伝子は、補正後に引き上げられて見える
次に、doubtサンプルで共通して発現上昇して見える遺伝子群（high@doubt）が、 各補正後の分布上でどの位置にあるのかを確認しました。 図5では、high@doubt遺伝子群を黒で示しています。
Gene Countデータでは、high@doubt遺伝子群はヒストグラムの左側に偏っています。 つまり、もともとのcount値としては低い領域に位置する遺伝子群です。
しかし、TMM-CPMでは、doubtサンプルのヒストグラム全体が右側に移動することで、 これらの遺伝子群が補正後に相対的に高く見えるようになっています。 その結果、high@doubt遺伝子群として、doubtサンプルで共通して発現上昇しているように見えます。
VSTやQuantile normalizationでは、TMM-CPMのように単純な横方向の移動としては見えません。 しかし、効果としては同様に、 もともと狭いダイナミックレンジの中にあった遺伝子群が補正後に引き上げられ、 doubtサンプルで相対的に高く見える状態が残っています。
low@doubt遺伝子は、補正の効果は見られるが完全ではない
次に、doubtサンプルで共通して発現低下して見える遺伝子群（low@doubt）について、 各データ上での分布を確認しました。 図6では、low@doubt遺伝子群を黒で示しています。
Gene Countデータを見ると、low@doubt遺伝子群は各サンプルのヒストグラムの右側、 つまり比較的高いcount値を持つ遺伝子群に位置しています。 しかし、doubtサンプルでは全体的なダイナミックレンジが狭いため、 これらの遺伝子のGene Count値が他のサンプルほど高くなっていません。 その結果、doubtサンプルでは相対的に発現が著しく低く見えます。
TMM-CPMでは、CPMスケールへの変換によってdoubtサンプルのヒストグラム全体が右側に移動し、 low@doubt遺伝子群の差は大幅に縮まりました。 ただし、差が完全になくなったわけではなく、 補正後もlow@doubt遺伝子群はやや低い値として残っています。
VSTやQuantile normalizationでは、TMM-CPMのように ヒストグラム全体が単純に横方向へ移動する関係は見えにくくなります。 しかし、これらの補正や変換でも、low@doubt遺伝子群の差は大幅に縮まっている一方で、 完全には補正しきれず、doubtサンプルでやや低く見える傾向が残っています。
青と緑の間にも、弱い分布の違いが残っている
doubtおよびclose to doubtのサンプルは、 全体のPCAやクラスタリングだけでも明確に他のサンプルと異なることが分かります。 一方で、high@doubtおよびlow@doubtの遺伝子群に注目すると、 青(4-mixed node)と緑(3-healthy node)のサンプル間にも、doubtサンプルほど大きくはないものの、 わずかな分布の違いが見られます。
このことは、サンプル構造が単純にdoubtサンプルだけで説明できるものではなく、 より弱い段階的な違いとして、残りのサンプル間にも存在している可能性を示しています。 そのため、doubtサンプルを除外した後も、 残りの青と緑のサンプル構造を確認する必要があります。
正規化や変換だけでは解決できない場合がある
ここまでの結果から、doubtサンプルに関連する発現パターンは、 raw Gene Countの段階でそのまま存在していた単純な発現差ではなく、 もともとのGene Count分布の歪みやダイナミックレンジの狭さが、 正規化・変換後の相対的な発現上昇または低下として現れたものと考えられます。
特に、ヒストグラムの形が大きく異なる、つまりサンプル間でダイナミックレンジが大きく異なる場合、 その影響は非常に厄介です。 TMM-CPM、DESeq2 VST、quantile normalizationのように異なる方法を試しても、 high@doubtおよびlow@doubtに対応する遺伝子群のパターンは残りました。 このことから、元データの分布の歪みは、 補正アルゴリズムを変えるだけで簡単に解決できる問題ではないことが分かります。
このように、生データの分布や発現プロファイルの形が大きく歪んでいる場合、 正規化方法や変換方法を変えるだけでは、 解析に適した状態に戻せないことはよくあります。 そのようなサンプルを無理に補正して解析に含めると、 DEG抽出やクラスタリングの結果を誤った方向に強調する可能性があります。
そのため、可視化によって明らかに他のサンプルとは異なるプロファイルを示すサンプルについては、 無理に補正して解析に組み込むよりも、 その理由をヒストグラム、PCA、クラスタリングなどの結果に基づいて明確に記録したうえで、 除外して解析を進める方が安全かもしれません。
歪みの強いサンプルを除外して、残りの構造を確認する
そこで次に、明らかに歪みの強いプロファイルを示す赤および黄色のサンプルを除外し、 残りの青および緑のサンプルだけを用いて再解析を行いました。
赤および黄色のサンプル群は、PCA、クラスタリング、ヒストグラムのいずれでも、 他のサンプルとは異なる発現プロファイルを示していました。 そのため、これらのサンプルを除外して再解析を行うことは、 可視化結果に基づく判断として妥当だと考えられます。
ただし、実験ノートや詳細なバッチ情報にアクセスできないため、 これらのサンプル群が、実際に不良サンプルであったのか、 特定の実験条件や技術的バッチに対応していたのかは判断できません。
ここでは、この除外後のデータを用いて、残りのサンプル構造を確認しました。
仮定したbatch情報でComBat補正を行うと、cluster間の差は弱まった
次に、明らかに歪みの強いサンプルを除いた後のデータについて、 残った緑と青のラベル情報を仮のbatch変数として扱い、ComBat補正を試しました。
ここで用いたラベル情報は、実験ノートやシーケンスrunなどに基づく 明示的なバッチ情報ではありません。 そのため、この解析は実際のバッチ効果を断定して補正するものではなく、 可視化上で見られた構造をbatchとして扱った場合に、 データがどのように変化するかを確認するための仮定的な検討です。
今回は、DESeq2 VST変換後の値に対してComBat補正を行いました。 VST値およびComBat補正後の値は、Subio Platform側で追加のlog2変換を行わず、 そのままクラスタリングに使用しました。
ComBat補正後のクラスタリングでは、 補正前に見られていたcluster間の大きな差はかなり弱まりました。 これは、今回ComBat補正を試した意図に沿う結果と言えます。
一方で、この結果だけから、 そのclusterが実際の技術的バッチであったと結論づけることはできません。 クラスタリングで見えた構造をbatchとして指定すれば、 その構造が弱まるのはある意味で当然です。 実際の研究データで同様の補正を行う場合には、 その構造を取り除いてよい根拠があるかどうかを必ず確認する必要があります。
DEGリストは補正方法によって変化する
ここまでの結果から分かるように、DEGリストはデータ構造や正規化・補正方法に依存して変化します。 逆に言えば、データの特徴や、正規化・補正に使った手法について正しく理解しないまま DEG解析を行うことは危険です。 DEGリストを解釈する前に、 その結果がどのようなデータ構造と前処理の上に成り立っているのかを確認する必要があります。
補正はデータを「正しくする」魔法ではない
今回の解析では、TMM-CPM、DESeq2 VST、そしてquantile normalization後のデータを比較しても、 doubtサンプルに関連する構造は残りました。 また、ヒストグラム上でも、high@doubtおよびlow@doubt遺伝子群の特徴的な位置関係は、 補正後も完全には消えませんでした。
一方で、明らかに歪みの強いサンプルを除いた後、 残った緑と青のラベル情報を仮のbatchとしてComBat補正を行うと、 補正前に見られていた緑と青のサンプル間の差は大きく弱まりました。 この結果は、補正方法によってデータの見え方やDEGリストが変化しうることを示しています。
しかし、それは単に「どの補正が正しいか」を選ぶ問題ではありません。 補正を行うには、何を補正対象とみなすのか、 それを取り除いてよい根拠があるのかを考える必要があります。
実験条件やバッチ情報にアクセスできない場合、 補正後の結果を正解として扱うのではなく、 複数の解析条件で安定して見える傾向と、 条件によって変わる部分を分けて考える必要があります。
正規化や補正は、データを「正しくする」ための魔法ではありません。 正規化後、変換後、補正後のデータを可視化し、 何が変わり、何が残ったのかを確認することが重要です。
歪んだデータを、後からアルゴリズムだけで「なかったこと」にすることはできません。 ComBatのような補正を使えば、指定した構造を弱めることはできます。 しかし、それは必ずしもデータが本来あるべき状態に戻ったことを意味するわけではありません。
最も重要なのは、解析後の補正に頼りきることではなく、 実験計画、サンプル調製、ライブラリ作製、シーケンス条件などの段階で、 できるだけ不要な揺らぎを抑えることです。 解析は、そのうえで得られたデータの状態を確認し、 補正してよい差と、補正してはいけない差を見極める作業です。